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Fourier  coefficients  of  target  signal  for  stationary 
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Fourier  coefficients  of  target  signal  for  non¬ 
stationary  target  at  frequency  oj^ 

Piecewise  linear  time  delay  parameters 
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Piecewise  linear  time  delay  trajectory  slope  parameters 

Smallest  time  interval  used  for  practical  correlation 
measurement 

Incremental  bandwidth  of  Fourier  coefficients 
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EFFICIENT  RECURSIVE  BATCH  TIME  DELAY 
DIFFERENCE  ESTIMATION  IN  THE  PRESENCE 
OF  TARGET  MOTION 


CHAPTER  1 
INTRODUCTION 

In  this  chapter  we  shall  define  the  problem  addressed  in  this  study 
along  with  our  associated  technical  objectives.  We  shall  give  back¬ 
ground  information  on  our  chosen  problem,  including  a  synopsis  of 
relevant  previous  works.  We  shall  state  clearly  the  organization  of  the 
material  presented  in  this  study.  Finally,  we  shall  make  a  statement  of 
contribution  along  with  an  executive  summary  of  our  major  results. 

1.1  BACKGROUND 

We  are  concerned  with  the  problem  of  estimating  time  delay  differ¬ 
ence  parameters  between  time  histories  of  broad  spectrum  energy 
observed  at  two  spacially  separated  sensors.  The  observations  at  each 
sensor  consist  of  measurements  made  on  a  combination  of  time-varying 
fields  of  coherent  broad  spectrum  energy  radiating  from  a  single  moving 
target  and  incoherent  local  broad  spectrum  noise  energy.  Each  measure¬ 
ment  consists  of  T  second  observations  of  the  received  time  varying 
energy  field. 
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There  are  numerous  estimators  and  associated  algorithms  designed 
to  estimate  the  time  delay  difference  between  two  sensors  due  to  a  sta¬ 
tionary  target.  However,  we  are  primarily  concerned  with  the  more 
difficult  problem  of  estimating  time  delay  difference  parameters  for 
significant  target  motion  at  low  Signal-power  to  Noise-power  Ratio 
(SNR).  In  the  real  world,  targets  of  interest  are  generally  in  motion 
and  initial  detection  of  targets  is  generally  made  at  low  SNR. 

The  above  problem  has  two  significant  phases.  The  first,  and  by 
definition  the  most  critical,  is  the  initialization  of  the  estimation 
process  on  target  until  a  steady  state  condition  is  achieved.  The 
second  is  the  continuing  steady  state  recursive  estimation  of  the 
target  time  delay  difference  parameters  for  the  duration  of  the 
problem.  In  the  past,  practical  estimation  algorithms  have  placed  most 
emphasis  on  the  later  steady  state  estimation  process  with  nominal 
concern  for  the  initialization  problem. 

The  traditional  approach  to  the  above  estimation  procedure  is  to 
use  recursive  adaptive  (digital)  filtering  techniques.  These  filter 
estimators  combine  the  current  measurement  with  the  previous  cycle 
filter  states  to  obtain  the  current  estimate.  Since  the  previous  filter 
states  are  a  function  of  all  previous  measurements,  these  filters  are 
called  Infinite  Impulse  Response  (HR)  filters.  The  detail  design  of 
these  filter  algorithms  represents  a  compromise  between  parameter 
estimate  variance,  dynamic  parameter  following  capability,  and  tran¬ 
sient  characteristics. 
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Our  approach  to  the  above  problem  shall  be  to  simultaneously 
operate  on  a  finite  block  of  measured  data  represented  by  M  consecutive 
T-second  observations.  This  type  of  estimator  is  called  a  batch  esti¬ 
mator.  If  the  finite  block  of  M  measurements  is  sequentially  updated 
with  new  measurements  every  T  seconds  while  discarding  the  oldest 
measurement,  the  estimator  is  called  a  recursive  batch  estimator.  In 
the  parallel  field  of  digital  filtering,  the  equivalent  filter  type  is 
called  a  Finite  Impulse  Response  (FIR)  filter. 

We  have  chosen  a  batch  estimator  for  several  reasons.  First,  main¬ 
taining  a  history  of  all  input  measurements  up  to  the  block  size  will 
yield  the  potential  for  good  initialization  algorithms.  Second,  a 
maintained  history  of  data  measurements  and  their  associated  time  delay 
difference  estimates  lend  themselves  naturally  to  a  priori  information 
assumptions  on  the  target  time  delay  difference  trajectory;  and  finally 
individual  measurements  in  the  measurement  block  can  respond  quickly  to 
target  dynamics  while  maintaining  the  total  observation  time  of  the 
entire  block. 

The  advantages  of  a  recursive  batch  estimator  are  two:  first,  a 
single  measurement  update  results  in  new  time  delay  difference  param¬ 
eter  updates;  and  second,  time  delay  difference  parameter  estimates  are 
based  solely  on  local  measurements. 
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Our  objective  is  to  investigate  and  design  optimal  recursive  batch 
estimators  for  the  two-sensor  non-stationary  target  time  delay  estima¬ 
tion  problem.  Our  estimation  algorithm  should  have  good  initialization 


properties  and  be  capable  of  estimating  significant  target  motion  at 
low  SNR. 

The  practical  problem  which  motivates  our  study  is  from  the  field 
of  underwater  acoustics.  We  want  to  estimate  the  bearing  and  bearing 
rate  (equivalent  to  time  delay  difference  and  time  delay  difference 
rate)  from  passive  spacially  separated  measurements  of  a  broad  spectrum 
noise  field  due  to  a  moving  target.  However,  the  recursive  batch  tech¬ 
niques  we  shall  examine  are  applicable  to  any  estimate  of  time  delay 
difference  due  to  measurements  of  broad  spectrum  signal  energy  from 
spacially  separated  point  sensors.  .In  addition,  the  recursive  batch 
algorithm  development  will  have  general  application  to  similarly  based 
algorithms. 

Practical  problems  require  practical  realizable  solutions.  There¬ 
fore,  our  work  will  strive’  to  apply  practical  considerations  to  theo¬ 
retical  work  at  every  opportunity.  In  addition,  a  major  objective  of 
this  study  is  to  develop,  analyze,  and  test  the  ground  work  for  candi¬ 
date  estimation  algorithms  based  on  our  theoretical  results. 

1.2  TECHNICAL  OBJECTIVES 

Major  technical  objectives  of  this  study  will  involve  three 
related  areas:  (1)  theoretical  development,  (2)  numerical  study,  and 
(3)  validation  via  computer  simulation. 


Our  principal  theoretical  development  will  be  to  derive  a  Maximum 
A  Posteriori  (MAP)  estimator  for  the  non-stationary  time  delay  estima¬ 
tion  problem.  A  key  step  in  this  process  is  to  develop  suitable  models 
for  the  non-stationary  time  delay  process.  From  our  derived  models  a 
Maximum  Likelihood  (ML)  estimator  and  a  priori  information  function 
will  be  derived  to  complete  our  MAP  estimator.  The  MAP  estimator  will 
be  extended  to  estimate  time  delay  difference  parameter  trajectories 
over  a  block  of  input  data.  The  time  delay  difference  parameters  to  be 
estimated  include  the  following: 

1.  Time  delay  difference  trajectory 

2.  Time  delay  difference  rate  trajectory 

3.  Time  delay  difference  acceleration  trajectory. 

Other  related  parameters  to  be  estimated  are: 

1.  Average  received  SNR 

2.  A  spectral  shape  parameter. 

Our  principal  numerical  effort  will  be  to  develop  an  efficient 
practical  algorithm  to  implement  our  MAP  estimators  over  a  block  of 
input  data  measurements.  The  practical  algorithms  will  take  on  two 
distinct  forms.  The  first  form  will  be  algorithms  using  discrete 
dynamic  programming  techniques.  The  second  form  will  be  algorithms 
using  gradient  search  techniques.  In  both  forms  our  algorithms  will 
stress  practical  application. 
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Our  first  experimental  objective  will  be  to  determine  the  perform¬ 
ance  characteristics  of  our  major  algorithms  as  a  function  of  algorith¬ 
mic  parameter  selection.  This  work  will  yield  insight  into  intelligent 
algorithmic  parameter  selection. 

Our  second  experimental  objective  is  to  study  relative  algorithm 
performance  and  performance  against  standard  algorithms  and  theoretical 
bounds. 

1.3  PREVIOUS  WORK 

The  problem  of  optimum  time  delay  estimation  under  the  assumption 
of  a  stationary  single  target  has  been  thoroughly  studied  by  many  inves¬ 
tigators.  Their  work  is  adequately  summarized  in  a  paper  by  C.  H.  Knapp 
and  G.  C.  Carter*  on  the  generalized  correlation  methods  for  estimation 
of  time  delay.  This  work  is  of  basic  interest  to  us  since  our  general 
approach  to  time  delay  estimation  is  also  using  correlation  techniques. 
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The  problem  of  optimum  time  delay  estimation  under  the  assumption 
of  a  single  non-stationary  target  has  received  less  attention. 

However,  some  papers  have  been  written  in  this  area.  Notable  among 
these  is  a  paper  by  Knapp  and  Carter  on  extending  the  generalized  cross 
correlation  techniques  to  the  case  of  a  moving  target.  Their  concept  is 
to  time  compress  or  expand  one  or  both  of  the  two  received  spacially 
separated  waveforms  to  compensate  for  target  motion  and  then  perform  a 
standard  generalized  cross-correl ation  technique  for  estimating  the 
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time  delay.  Although  the  Knapp  and  Carter  work  leads  to  a  signifi¬ 
cantly  different  implementation  strategy  than  our  work,  we  shall  make 

extensive  use  of  their  general  approach  to  calculate  our  ML  estimator. 

3 

In  addition,  Y.  T.  Chan,  J.  M.  Riley,  and  J.  B.  Plant  have  done  work 

where  the  non-stationary  time  delay  process  is  viewed  as  a  FIR  filter, 

4 

and  Friedlander  has  done  work  in  combined  spectral  and  time  delay  esti 
mation  using  an  Auto  Recursive  Moving  Average  (ARMA)  model.  The  time 
delay  estimation  portion  involves  a  moving  average  or  FIR  filter 
representation  of  the  unknown  non-stationary  time  delay.  We  have  made 
less  use  of  the  above  work  in  the  development  of  our  estimators. 
However,  the  FIR  model  for  the  non-stationary  time  delay  process  has  a 
parallel  in  our  development  of  a  'MISMATCH'  function. 

5 

L.  C.  Ng  in  his  PhO  dissertation  addresses  the  problem  of  multi¬ 
target,  multisensor  tracking  of  a  non-stationary  target.  The  batch 
processing  approach  and  general  time  delay  trajectory  modeling  used  by 
Ng  is  similar  to  our  own  development.  We  shall  expand  on  the  work  of 
Ng  for  the  two-sensor  one-moving-target  scenario.  Our  most  significant 
development  in  this  area  shall  be  the  improved  derivation  of  the 
'MISMATCH'  function  and  its  redefinition  for  use  in  a  practical 
estimator. 

The  general  approach  for  practical  time  delay  estimation  involves 

maximizing  a  MAP  function  over  a  block  of  data.  An  important  practical 

implementation  of  the  above  approach  was  inspired  by  the  dynamic  pro- 
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gramming  technique  employed  by  J.  J.  Wolcin  ’  ’ 


and  refined  to  a 


recursive  batch  process  by  R.  A.  LaTourette,  S.  G.  Greineder,  and 
*  9 

J.  J.  Wolcin. 

1.4  REPORT  ORGANIZATION 

The  organization  of  this  report  basically  follows  our  technical 
objectives.  Chapter  1  defines  our  problem  and  summarizes  our  contribu¬ 
tion  and  results.  Chapter  2  discusses  our  assumed  models.  Chapter  3 
develops  the  theoretical  MAP  estimator.  Chapter  4  describes  our  prac¬ 
tical  algorithm  development.  Chapter  5  gives  the  results  of  our  empir¬ 
ical  algorithm  parameter  investigation.  Chapter  6  presents  the  results 
of  our  relative  algorithm  performance  studies.  Finally,  Chapter  7  sum¬ 
marizes  our  work  and  suggests  future  work. 

1.5  STATEMENT  OF  CONTRIBUTION  ANO  EXECUTIVE  SUMMARY 

We  have  derived  an  approximate  MAP  estimator  for  the  time  delay 
difference  parameters  of  a  moving  acoustical  target  as  seen  by  two 
sensors.  The  target  is  assumed  to  radiate  flat  broad  spectrum  Gaussian 
noise.  The  ML  portion  of  the  MAP  estimator  was  determined  to  be  the  set 
of  parameters  which  maximized  the  product  of  the  measured  normalized 
cross-correlation  lag  points  between  point  sensors  and  the  normalized 
expected  cross-correlation  lag  points  (mismatch  function). 

For  the  case  of  linear  time  delay  motion,  we  have  calculated  the 
correct  analytical  form  of  the  mismatch  function  and  the  appropriate 
normalization  factor.  To  the  best  of  our  knowledge,  this  is  the  first 


correct  derivation  of  the  mismatch  function  to  appear  in  the 
literature. 

The  derivation  of  the  mismatch  function  was  expanded  to  include 
the  more  general  case  of  a  linear  target  power  spectrum.  From  this  work 
we  were  able  to  develop  an  estimator  for  a  parameter  related  to  the 
slope  of  the  unknown  target  power  spectrum. 

For  our  ML  estimator  we  have  assumed  that  the  received  target  time 
delay  trajectory  is  linear  over  our  measurement  interval  and  continuous 
(piecewise  linear)  between  measurements.  However,  over  a  block  of 
M  contiguous  measurements  we  have  assumed  that  the  received  target  time 
delay  difference  trajectory  is  a  quadratic  function  of  time.  From  the 
later  assumption  we  have  derived  an  a  priori  information  penalty  func¬ 
tion  based  upon  the  sum  of  the  squared  residuals  between  the  candidate 
ML  time  delay  difference  estimates  and  their  best  quadratic  least  mean 
square  fit  (LMSF).  The  role  of  the  LMSF  a  priori  information  penalty 
function  is  to  enhance  the  local  ML  estimate  with  information  obtained 
over  the  whole  block  of  M  measurements.  We  noted  that  using  the  LMSF 
a  priori  information  penalty  function  simultaneously  with  the  ML  esti¬ 
mate  was  different  and  superior  to  post  refining  the  ML  estimates  with  a 
quadratic  LMSF.  The  quadratic  LMSF  penalty  function  also  yields  a 
procedure  for  obtaining  refined  estimates  of  time  delay  difference, 
time  delay  difference  rate,  and  time  delay  difference  acceleration. 

The  optimum  relative  weight  K  between  the  ML  estimate  and  the  a  priori 
information  was  determined  to  have  a  complex  functional  dependency 
including  algorithm  implementation. 


A  major  emphasis  of. our  study  was  to  derive  practical  recursive 
block  estimation  algorithms  based  on  our  theoretical  work.  Therefore, 
our  implementation  problem  can  be  defined  as  follows:  we  wish  to 
determine  the  time  delay  difference  parameters  which  maximizes  our 
derived  MAP  function  over  a  block  of  M  consecutive  measurements.  We 
further  require  that  the  estimates  of  time  delay  difference  parameters 
be  recursively  updated  as  the  oldest  data  measurement  is  removed  from 
the  data  block  and  a  new  measurement  is  added  to  the  data  block. 

For  the  above  objective,  we  have  derived  an  efficient  recursive 
dynamic  programming  algorithm  and  a  recursive  gradient  search 
algorithm. 

The  basic  dynamic  programming  algorithm  is  based  upon  a  similar 
algorithm  developed  by  Wolcin^’^’^  for  frequency  estimation.  We 
claim  credit  for  introducing  the  recursive  modification  to  the  algo¬ 
rithm,  modifying  the  algorithm  for  the  time  domain  application,  and 
increasing  algorithm  efficiency. 

The  recursive  modification  for  the  dynamic  programming  algorithm 
involved  recursively  updating  the  sum  of  square  residuals  of  a  quad¬ 
ratic  LMSF.  For  this  purpose  we  employed  an  algorithm  developed  by  Ng^ 
for  efficiently  updating  recursive  least  squares  parameters.  In  study¬ 
ing  Ng's  algorithm,  we  discovered  a  numerical  instability  and  imple¬ 
mented  a  numerical  technique  to  eliminate  this  instability. 


The  gradient  search  algorithm  is  based  upon  calculating  the 
gradient  vector  of  the  MAP  function  over  the  block  of  input  measurements 
with  respect  to  the  unknown  time  delay  difference  initial  offset  and 
rates.  We  have  analytically  calculated  the  required  gradient  vector. 

The  algorithm  uses  the  previous  update  parameter  estimates  to  determine 
excellent  initial  parameter  estimates  for  the  current  measurement 
update.  To  the  best  of  our  knowledge,  this  is  the  first  numerical 
search  technique  to  be  applied  to  this  type  of  problem. 

For  both  the  dynamic  programming  and  gradient  search  algorithms  we 
have  recognized  and  have  placed  special  emphasis  on  improving  initial 
acquisition  properties. 

The  following  is  a  top  level  list  of  major  results  obtained  from 
simulation  experiments.  A  more  detailed  list  can  be  found  in  Chapter  7. 

1.  The  gradient  search  algorithm  was  found  to  have  superior 

steady  state  performance  characteristics  with  respect  to  other 
related  recursive  batch  algorithms  and  the  Cramer-Rao  Lower 
Bound  (CRLB)  of  a  coherent  integrator  with  equivalent  averag¬ 
ing  time  but  no  a  priori  information.  The  gradient  search 
algorithm  obtains  its  advantage  over  the  CRLB  by  the  addition 
of  a  priori  information. 


2.  The  gradient  search  algorithm  yields  meaningful  estimates  of 
both  time  delay  difference  and  time  delay  difference  rate. 

The  time  delay  difference  acceleration  estimate  was  not  mean¬ 
ingful. 

3.  The  gradient  search  algorithm  only  yields  a  valid  target  power 
spectral  shape  parameter  estimate  for  SNR  level  greater  than 
or  equal  to  -10  dB. 

4.  Moderate  target  power  spectral  modeling  errors  have  little 
affect  on  the  performance  of  the  gradient  search  algorithm. 

5.  Optimal  parameter  selection  can  be  determined  for  the  gradient 
search  algorithm.  However,  the  optimum  parameters  are  often 
dependent  on  the  target  motion  characteristics. 

6.  The  dynamic  programming  algorithm  has  excellent  initialization 
properties.  With  moderate  target  dynamics,  the  dynamic  program 
algorithm  can  successfully  initialize  within  a  few  dB  of 
steady  state  performance  thresholds.  The  gradient  search 
algorithm  initialization  performance  is  slightly  degraded  from 
the  dynamic  programming  algorithm.  The  gradient  search 
algorithm  initialization  performance  is  still  very  good,  but 
there  is  potential  for  improvement. 


CHAPTER  2 


PROBLEM  FORMULATION 


2.1  INTRODUCTION 

In  this  section  we  shall  discuss  the  assumed  observation  scenario, 
mathematically  describe  the  resulting  input  sensor  measurements,  and 
formulate  the  general  batch  Maximum  A  Posteriori  (MAP)  time  delay 
difference  estimation  problem.  Since  a  major  objective  of  this  study  is 
to  realize  practical  estimators,  practical  considerations  shall  be 
noted  and  incorporated  throughout. 

2.2  DESCRIPTION  OF  THE  OBSERVATION  SCENARIO 

Let  us  assume  a  sensor  array  which  can  be  represented  as  two 
sensors  in  space  separated  by  a  distance  Dp  (Dipole  Spacing).  Let  us 
further  assume  that  a  distant  target  has  relative  motion  with  respect  to 
the  received  sensor  array.  The  target  is  at  a  sufficient  range  from  the 
sensor  array  that  the  received  wavefronts  can  be  considered  plane 
waves.  Figure  2-1  depicts  the  scenario  below.  The  following  equation 
states  the  condition  for  which  our  plane  wave  assumption  is  valid. 

L2 

RANGE  >  >  ^2.  (Steinberg) 11  (2-1) 

Lg  =  Array  or  aperture  length 

X  =  Wave  length  of  tne  highest  frequency. 
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Figure  2-1.  Observation  Scenario 


In  practice,  if  the  target  is  a  large  object  such  as  a  ship,  or  to 
a  lesser  extent  an  aircraft,  the  motion  of  the  target  will  be  continuous 
and  slowly  changing  over  an  extended  period  of  time.  We  can  assume  that 
the  resulting  relative  bearing  trajectory  and  the  corresponding  time 
delay  trajectories  can  be  adequately  modeled  by  a  low-order  polynomial 
over  the  time  interval  of  interest. 


■*» 


2.3  DESCRIPTION  OF  INPUT  SENSOR  MEASUREMENTS 


From  Figure  2-1  the  received  signal  at  both  sensor  points  can  be 
modeled  as  shown  below. 


x  1  ( t )  =  s(t  +  o^t);  +  n1  i  t ) 

X9(t)  =  s ( t  +  ;,(t })  +  n9(t ) 


where  s(t)  is  the  signal  waveform  and  r^(t)  and  n2(t)  are  the  noise 
fields  at  sensor  1  and  sensor  2,  respectively. 


We  shall  assume  that  s(t),  n ^ ( t ) ,  n2(t)  are  zero  mean,  Gaussian  and 
mutually  uncorrelated  random  processes  with  known  band-limited  flat* 
power  spectrums;  i.e., 


Gs(f)  = 


2 ( f 2  -  f1)t 


fL  <  f  <  f2  or  -f2  <  f  <  -f1 


0, 


otherwise 


(2-3) 


G„  (f)  3  G„  (f)  - 


N 


'1 


2(f,  -  M’  T11T-T2 

otherwise 


2  '  'V 

0, 


fi  <  f  <  fo  or  -f2  <  f  <  -fx 


(2-4) 


where 

As  =  Total  signal  power 
N$  =  Total  noise  power 
f^  =  lower  bandpass  cutoff  frequency 
f^  =  upper  bandpass  cutoff  frequency. 


*  Although  we  assume  that  both  the  signal  and  noise  spectrums  are  flat, 
in  practice,  this  assumption  will  rarely  be  true.  The  normal  prac¬ 
tice  is  to  prewhiten  or  flatten  the  anticipated  noise  spectrum.  The 
nominal  assumption  is  that  the  signal  spectrum  will  have  the  same 
shape  as  the  noise  spectrum  and  thus  our  prefilter  will  be  optimal. 
Later  in  the  analysis  we  shall  loosen  the  requirement  that  the 
filtered  signal  spectrum  will  be  flat  and  develop  an  algorithm  to 
estimate  a  signal  spectrum  parameter.  However,  the  primary  develop¬ 
ment  shall  be  with  a  flat  signal  and  noise  spectrum  assumption. 
Therefore,  emphasis  must  be  placed  on  analyzing  performance  wit-i  non¬ 
flat  spectrums. 
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In  the  above  definition  there  are  two  implied  assumptions.  The 
first  is  that  the  nominal  receive  signal  power  is  identical  at  sensors  1 
and  2.  The  second  is  that  the  noise  processes  and  the  signal  process 
have  identical  bandpasses.  In  practice,  both  assumptions  are  not 
unreasonable.  In  general,  we  would  expect  identical  signal  power  at 
both  receive  sensors.  Also  for  most  practical  applications  the 
received  signals  are  processed  through  some  bandpass  filter  with  a 
smaller  bandwidth  than  either  the  signal  process  or  the  noise  process. 
Therefore,  the  cutoff  frequencies  of  the  input  bandpass  process  can  be 
considered  as  the  limiting  frequencies  f^  and  f2> 

<(>^(t)  and  <t>2 ( t )  represent  the  received  signal  delay  as  seen  by 
sensor  1  and  sensor  2,  respectively,  due  to  the  moving  target  as  a 
function  of  time. 


We  shall  assume  that  over  the  total  observation  interval  (r)  that 
and  <}>2(t)  can  be  adequately  modeled  as  a  low  order  polynomial  or 
specifically  as  a  quadratic  function;  i.e.. 


♦x(t)  =  d*  t2  +  dj  t  +  dj 

<t>2(t)  =  d2  t2  +  d2  t  +  d2 


0  <  t  <  T 


0  <  t  <  T 


(2-5) 


where  the  d's  are  coefficients  of  a  polynomial  representing  time  delay 
as  a  function  of  time. 
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We  shall  further  assume  that  if  the  total  observation  interval  r  is 
subdivided  into  M  equal  subdivisions  each  of  length  T,  that  <t>^(t)  and 
cf>2 ( t )  can  also  be  adequately  modeled  as  a  piecewise  linear  function 
•defined  below  for  4>^(t). 

4>i(t)  =  <t>i_1(tk_1)  +  AsJ(t  -  tk_1),  tk-1  <  t  <  tk, 

k  =  1,  2,  3,  ...  M  (2-6) 
4*5(0)  =  dj  (2-6A) 

where  now  we  define 

k  =  subinterval  number 

dj  =  initial  or  reference  time  delay  at  sensor  1 

AsJ  =  (dj  -  dk_x )/T  (2-6B) 

d£  =  the  time  delay  at  sensor  1  at  time  tk. 

Note  that  for  this  model  the  d's  now  represent  actual  time  delays 
at  given  time  instances. 

2.4  STATEMENT  OF  PROBLEM 

Our  objective  is  to  investigate  optimal  and  suboptimal  estimators 
for  the  time  delay  difference  between  two-point  sensors  for  a  single 
moving  target  as  a  function  of  time.  The  estimators  we  will  examine 
shall  operate  on  M  input  measurements,  each  of  T  second  duration,  from 
the  two-point  sensor  array.  The  measurements  will  oe  f^om  the  current 
time  back  to  M  •  T  or  x  seconds  in  the  past.  This  type  of  process  is 
known  as  a  batch  process  and  our  estimators  will  be  known  as  batch 


\ 
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estimators.  The  batch  process  has  the  advantage  of  receiving  updates, 
processing  data,  and  making  decisions  at  a  T  second  rate  while  maintain¬ 
ing  a  total  observation  time  of  M  •  T  or  x  seconds.  A  primary  objective 
for  our  estimators  will  be  to  enhance  the  above  advantage. 

In  general,  a  batch  estimator  yields  estimates  for  any  point 

5 

within  the  processed  block  of  data.  The  most  accurate  estimate  is  at 
an  interior  point  of  the  input  block,  while  the  least  accurate  estimate 
is  at  the  end.  For  our  study  we  shall  be  interested  in  the  first  or  most 
recent  estimate  and  the  middle  estimate. 

A  second  requirement  for  our  approach  is  that  they  are  recursive 
estimators  (i.e.,  for  every  new  sensor  measurement  (T  seconds),  there 
will  be  new  estimator  outputs).  This  requires  that  our  input  block  of  M 
sensor  measurements  acts  as  a  sliding  window.  When  a  new  sensor 
measurement  is  received,  it  is  added  to  the  block  of  input  data  while 
the  oldest  sensor  measurement  is  dropped  from  the  block  of  input  data. 

New  estimates  will  be  calculated  after  each  update.  For  the  case  of  the 
middle  time  frame  estimate,  this  process  will  result  in  a  time  history 
of  estimates  each  separated  by  T  seconds  in  time  and  the  most  recent 
estimate  delayed  by  one  half  the  block  length.  The  most  recent  time 
frame  estimate  will  also  result  in  a  time  history  of  estimates  each 
separated  by  T  seconds,  but  the  most  recent  estimate  will  be  nominally 
at  current  time.  Obviously,  with  the  appropriate  time  shift  of  one  half 
the  data  block  length,  these  two  time  histories  will  coincide,  we  would 
expect  that  the  variance  of  the  middle  t’ne  frame  estimates  would  :e 


less  than  the  most  recent  time  frame  estimates. 


As  mentioned  in  Chapter  1,  we  will  approach  the  above  problem 
using  MAP  techniques.  The  MAP  approach  was  choosen  to  take  advantage  of 
a  priori  knowledge  of  the  characteristics  of  the  target  time  delay 
difference  trajectory.  This  analysis  will  yield  the  form  of  the  input 
sensor  measurements  and  optimum  or  near  optimum  processing  strategy. 
These  results  will  be  directed  towards  practical  considerations  to 
yield  a  basis  for  practical  estimators. 


CHAPTER  3 


THEORETICAL  DEVELOPMENT  OF  MAP  ESTIMATORS 

3.1  INTRODUCTION 

In  this  chapter  we  shall  examine  the  two-sensor,  single  moving 
target  time  delay  difference  estimation  problem  outlined  in  Chapter  2. 

Our  objective  is  to  determine  a  form  for  the  input  measurements  and 
derive  an  optimum  or  near  optimum  processing  strategy.  A  key  concept  is 
that  we  are  looking  at  two  problems  in  one.  On  one  hand  we  would  like  to 
view  our  problem  as  making  a  single  measurement  over  a  large  interval  of 
time  t.  For  this  problem  we  assume  that  the  unknown  time  delay  differ¬ 
ence  trajectory  is  quadratic  in  time.  This  approach  has  the  advantage 
of  long  integration  time  of  the  observed  measurements.  We  shall  refer 
to  the  above  as  Problem  1.  A  second  view  is  that  our  measurements  con¬ 
sist  of  M  intervals  of  duration  T  seconds,  where  M  •  T  =  x.  The  total 
time  delay  difference  trajectory  is  assumed  to  be  piecewise  linear  with 
each  of  the  M  subintervals  being  linear.  This  approach  has  the  advan¬ 
tage  of  being  more  adaptive  (T  second  update  rate)  to  the  time  delay 
difference  trajectory  dynamics  and  has  a  simpler  model  (linear)  than 
Problem  1.  We  shall  call  this  Problem  2.  Ideally,  we  would  like  to 
simultaneously  gain  the  advantages  of  long  observation  time  of  the 
measurement  data  of  Problem  I  and  the  adaptability  of  Problem  2. 
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3.2  THE  MAP  EQUATION 

Our  general  approach  shall  be  to  apply  Maximum  A  Posteriori  (MAP) 
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estimation  techniques  as  outlined  by  Van  Trees.  The  MAP  estimate  is 
defined  as  the  set  of  parameters  a(O)  which  maximizes  the  a  posteriori 
density  function  Pa/Q(A/0)  of  a  random  parameter  vector  A  conditioned 
on  observations  or  measurements  0. 

The  standard  procedure  is  to  rewrite  pa/o(A/0)  in  terms  of  the 
measurements  0  conditioned  on  the  unknown  parameter  vector  A  (maximum 
likelihood  estimate)  and  a  priori  knowledge.  By  Bayes  rule,  we  have  the 
fol lowing. 


Po/,( 0/A)  P.(i) 
fa/o(A/0)  *  0/3  Po(0)a 


(3-1) 


Taking  the  logarithm  of  Equation  3-1  which  will  not  alter  the 
location  of  the  maximum  yields  the. fol lowing  result: 


In  pa/o(A/0)  =  ^  PQ/a(0/A)  +  In  Pg(A)  -  In  PQ(0) 


(3-2) 


Since  the  measurement  space  0  does  not  depend  on  the  parameter 
vector  A,  the  final  MAP  equation  becomes  the  following: 


MAX  [In  Po/a(0/A)  +  In  P#(A)] 
A 
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Equation  3-3  contains  two  components.  The  first  [In  PQ/a(0/A)] 
represents  a  Maximum  Likelihood  (ML)  estimate  of  the  parameter  vector  A 
from  the  observed  measured  data  0.  The  second  component  [In  P _ (A) ] 

—  a  *■— 

represents  a  priori  knowledge  of  the  unknown  parameter  vector  A. 

In  the  following  sections  our  goal  will  be  to  develop  expressions 
representing  the  ML  and  a  priori  portions  of  the  MAP  estimator. 

3.3  THE  MAXIMUM  LIKELIHOOD  ESTIMATOR  FOR  PROBLEM  1 

Our  objective  is  to  find  the  form  of  the  ML  estimator  for  the  time 
delay  difference  trajectory  between  two-point  sensors  for  a  single 
moving  target.  We  assume ‘that  we  observe  the  sensor  measurements  for  a 
period  of  time  t  and  that  the  received  signal  time  delay  as  seen  by  each 
sensor  can  be  represented  by  a  quadratic  function  in  time.  Therefore, 
one  can  describe  the  time  delay  trajectory  at  sensors  1  and  2,  respec- 
tively,  as  the  parameter  vectors  D  (d2>  d^,  dQ)  and  D  (d^,  dj,  dQ) 
which  represent  coefficients  to  a  quadratic  equation  in  time.  The  time 
delay  difference  trajectory  AD(Ad2>  Ad^,  Adg)  is  also  quadratic  in  time 


Equation  2-5  for  the  time  delay  difference  A<j>(AD,  t)  is  repeated 
below  for  convenience. 


A4>(AD,  t)  3  Ad2  t2  +  Adx  t  +  AdQ  .  (3-5) 

Our  approach  to  initially  estimate  the  form  of  the  ML  estimator  for 
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Problem  1  will  follow  the  lines  developed  by  Knapp  and  Carter. 


We  have  from  Chapter  2  and  above  the  following  models  for  the 
signal  received  at  sensor  1  and  sensor  2,  respectively. 


Xj_(t)  3  s (t  +  <fr1(D.1,  t))  +  n^t) 
x 2 ( t )  3  s (t  +  <j>2(02,  t) )  +  ri2(t) 
where  s(t)  and  n(t)  are  described  in  Chapter  2. 


(3-6) 


The  Knapp  and  Carter  approach  is  to  represent  the  observed  signal 
waveforms  [x^(t),  x2(t)]  as  a  set  of  Fourier  coefficients 

^f^-rP’  Cf^aJ-N+l^  cf^tV’  Cf2^-N^’  Cf2^-N+1^  *•'  cf2^N^‘ 
These  Fourier  coefficients  are  assumed  to  be  independent  zero  mean 

Gaussian  random  processes.  This  assumption  assumes  a  sufficient  Band¬ 
width  Time  (BT)  product  (i.e.,  BT  >  8).^  The  procedure  for  calculating 


the  ML  estimate  is  to  calculate  the  joint  conditional  Gaussian  density 
function  of  the  Fourier  coefficients  conditioned,  at  least  initially, 

-1  ~p 

on  the  unknown  parameter  vectors  D  ,  D  .  The  ML  estimate  is  the  vectors 


0  ,  of  which  maximize  the  log  of  the  conditional  Gaussian  density  func¬ 
tion.  However,  it  is  important  to  note  that  our  primary  interest  is  to 
estimate  the  unknown  time  delay  difference  parameter  vector  AD. 


Proceeding  along  the  above  line  of  reasoning,  we  can  represent  the 
kth  frequency  component  of  x^(t)  and  X2U)  as  the  following: 


(3-7) 


where  1^  !  k  Hoi  and  Ad  3  2ir/x  . 

We  can  now  represent  our  observation  vector  [x-^(t),  X2 ( t )  ]  over 
the  interval  t  by  a  vector  yj  of  their  Fourier  coefficients  given  below: 

IT  =  [cJ(w_N),  cJ(d_n+1)  ....  cJ(uq)  ...  cJ(uN)]  (3-8) 

where  C^u^)  =  [cf  (0^),  cf  (c^)]  . 


(3-9) 


We  shall  now  roughly  determine  the  form  of  the  covariance  matrix  of 
Equation  3-12.  This  calculation  will  allow  us  to  obtain  the  general 
form  of  a  practical  approximate  ML  estimator.  In  later  analysis  we 
shall  re-examine  this  calculation  in  more  accurate  detail  with  a  dif¬ 
ferent  theoretical  model.  This  new  analysis  will  yield  the  precise  form 
of  the  'mismatch  function'. 


Our  approach  to  calculate  the  covariance  matrix  (Equation  3-12) 
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shall  be  similar  to  that  adopted  by  Ng. 


From  Equations  3-6  and  3-7,  we  have  that 


:fi(u»k)  ■  \  f  [s(t,  ^(D1,  t))  +  nx ( t ) ] e  ^dt  .  (3-13) 


The  above  expression  can  be  rewritten  as  a  signal  component  a^(uik) 
and  a  noise  component  N^(wk). 


cf  (wk)  =  0^(0^)  +  N^u^)  . 


(3-14) 


(oj^ )  is  the  Fourier  coefficient  at  frequency  ^  =  k  io  for  the 
noise  process  as  seen  at  sensor  1. 


c^(wk)  represents  the  Fourier  coefficient  at  ^  =  k  for  the 
non-stationary  signal  process  as  seen  at  sensor  1. 


^  )  “  T  /  $  (  t 


~1  -J-kt 
id  ,  t)e  dt 


.3-15; 


'•  .*’**•/•  AA  A  .*•  _*.'.*.**»  „  •  *  . -*A"  \*  s' ■  *  \*  \*  •.*  * 


We  would  like  to  be  able  to  express  the  Fourier  coefficient 

for  the  non-stationary  target,  in  terms  of  a^(ojk),  the  Fourier  coeffi¬ 
cient  for  a  stationary  target,  which  we  have  assumed  knowledge  of. 

The  stationary  Fourier  coefficient  for  sensor  1  is  given  below. 


and  the  stationary  band-limited  signal  s(t)  can  be  written  in  terms  of 
the  Fourier  coefficients  as  the  following: 


juj0t 


S(t)  =  2  ajU^e  £ 
a*-N 


(3-17) 


rKl 


In  Equation  3-15,  make  the  substitution  y  *  t  +  <j>^(D  ,  t), 
yielding: 


(O,1 ,  t) 


“l(u,k}  =  T  f  ,  $(^ 

’'As .  (  D  .  T 1 


-JVi(D  ,y) 


-^r 

*  t  n 


1  +  0, (D  ,  y) 


(3-13) 


Substituting  Equation  3-17  into  Equation  3-13  and  exchanging  sum¬ 


mation  with  integration  yields: 

N 

0^(0^)  =  ^  1  “l^)  1 1 ^ k T» 
£=-N 


where 


1  +  y) 


(3-19) 


(3-19A) 


We  can  repeat  the  above  formulation  for  sensor  2  to  yield  the 
fol lowing : 

N 


<*2 ( wk )  ”  a2 ( )  *2^ ’  —  ‘ 


(3-20) 


£»-N 


where 


1  ,T^  j  [Y^-t^^^2’^)] 

"  T  /  „  2  ^ 

•/4>2(r,-r) 


_dx. 


1  +  $?(D2,  y) 


(3-20A) 


We  now  have  the  equations  to  symbolically  calculate  the  four  elements  of 
the  covariance  matrix  of  Equation  3-12. 


The  (1,1)  element  using  Equations  3-12,  3-9,  ana  3-14  can  be 


written  as  the  following: 


E[cf  (u)k)  cf^(ojk)]  =  E  [a^(ojk)  Q^lu^)]  +  EfN^u^)  N^a^)]  .(3-21) 


Note  that  due  to  our  uncorrelated  processes  assumption. 


Efa^u^)  N ^ ( 00^ ) ]  =  E [Nj(ajk)  0^(0^)]  =  0 


(3-22) 


Letting  represent  the  expectation  operation  and  substituting 

Equation  3-19  into  Equation  3-21  yields  the  following. 


E[cf  (wk)  c f  (u>k ) ]  =  (a^u)^)  a^u^)  1^  1^)  +  N^u^)  N^(u>k )  . 

1  1  l-N  (3-23) 


If  we  now  let  x  effectively  approach  infinity  while  maintaining  cuk 
constant,  we  can  rewrite  Equation  3-23  as  an  integral.  Note  the  follow 
ing  relationships. 


_  2irl 

‘  T 


.  2ttN 
JN  '  T 


-2ttN 
-N  "  x 


(3-24A) 


(3-243) 


(3-24C) 


x  ajtu.^)  =  GsU£) 


(3-240) 


L  i  m  x  N 1  (  jj.  )  N  i  ( o>,  )  -  G  ( u,  ) 


,3-2 4E; 


->  J* 


where  G  (u)  is  the  signal  power  spectrum  described  in  Chapter  2  and 
(co)  is  the  noise  power  spectrum  as  seen  at  sensor  1  (also  described 
in  Chapter  2) . 


With  the  above  in  mind.  Equation  3-23  can  be  written  as  the 
fol lowing. 


□ 

E[Cf  U),  c*  (cuk)]  =  ^  f  G  (ai)  Limdj 
11  JQ  T*» 


(oi.  ,  CO,  T,  D1) 


CO,  t,  D^Jdco  +  ^  (cok)  . 


(3-25) 


Certainly  if  Gs(co)  is  a  constant  as  we  initially  assumed  in 
Chapter  2,  we  can  rewrite  Equation  3-25  as  the  following.  Note  that  we 
make  the  above  assumption  for  convenience  at  this  point.  Later  we  shall 
relax  this  restriction. 


E[cf  K)  cfi(cok)]  =  Gs(cok)  Q^x,  D1)  +  — ^ 


(3-26) 


where 


Q]_(wk,  T,  D1) 


1  f  ~1  *  -1 

/  Cim(I^( tok ,  co,  x,  J3  )  I  -  ( *  co »  x ,  0,  )dco 

J  a  1 


(3-26A) 


One  can  view  the  function  Q(-)  as  a  distortion  of  the  autopower 
spectrum  for  a  stationary  target  due  to  target  motion.  Tne  distortion 
will  be  in  terms  of  spectral  spreading  or  contraction  ana  magnitude 


•<  *v 


distortion.  For  reasonable  target  motion,  we  would  expect  the  auto¬ 


power  spectrum  distortion  to  be  small. 


Two  points  are  important  to  note  at  this  time.  The  first  is  that, 

at  low  SNR  where  the  noise  power  will  dominate,  small  signal  power 

losses  due  to  motion  will  be  impossible  to  separate  from  other  causes  of 

signal  power  loss.  The  second  point  is  that  sensor  1  autopower  spectrum 

depends  only  on  parameters  D*.  Our  primary  interest  is  to  maximize  with 

-2  ~1 

respect  to  the  difference  parameters  AD  =  _D  -  _D  . 


We  can  use  identical  arguments  to  obtain  an  expression  for 
★ 

E(cf  (u)k),  Cf  (cuk ) ]  or  element  (2,2)  which  is  written  below. 
t[cf^(u)k),  cf^(wk)]  =  Gs(wk)  Q2(uk,  t,  D2)  +  G^(u;k)  • 


(3-27) 


Using  similar  arguments  as  above  and  noting  that  E  [N^(oj).  N2(uj) ]  =  0 
we  can  write  an  expression  for  E[cf  (u>k)  cf  (wk)]  or  element  (1,2) 
which  is  shown  below. 


u 

E[cf1K)  cf2K)]  =  h / 


G  (u>.  )  Lim(I,(u>,  u)k,  t,  D  )  • 


I2(u),  u>k,  t,  5  ))dcj  .  (3-2S) 


Again  for  convenience  assuming  that  G$(u:k)  is  a  constant,  we  have 


t[Cr  (U)J,  c*  (u),  )  ]  -  GIujJ  Q,  0U.,  t,  D2,  D1)  (3-29) 


AV'.'VV 


W  W  ?  V  ‘  -  ' 


where 


o 

2^\'  5.  >  R  )  4^  f  Lim(I^(u],  t,  u^,  t,  D_  ))duj  • 

-B  ^ 

(3-29A) 

Finally,  we  can  write  E  [cf  (u>k),  c*  (u>k ) ]  or  element  (2,1)  as  the 
fol lowing: 


E  tCf 2 (^k )  ’  cf1^a)k^  =  GsH^  ^  1 , 2  ^  wk  ’  T’  ^ ^ 


(3-30) 


Summarizing  Equations  (3-26),  (3-27),  (3-29),  and  (3-30),  we  have 
the  following  matrix: 


Z((jJk)  =  t  E[Cf(<Dk)  Cf((jjk)] 

GS(W<J  T  Q1(^k,  t,  D1)  +  Sn^(wk)  G$(u)k)  T  Gl,2^k’  T’  £2) 

_GsUk)  T  Qi^^k’  T’  ~2^  Gs^wk^  1  G2^k’  T’  +  Gn  ^k^ 

(3-31) 


Z(^). 


Looking  ahead,  we  can  write  an  expression  for  the  inverse  of  matrix 
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2(“v>  ■  TzT^TT  • 

Gs(ujk)  t  Q2(u)k,  t,  £2)  +  Gn^(u>k)  -G^u^)  T  ^l,2^k’  T’ 

-65(00^)  t  T»  1^’  Ss^k^  T  Ql^“k»  T*  £*)  +  Gn  (t^) 

(3-32) 

where  the  determinant  |z(wk)|  can  be  written  as: 

\Z(^)\  =  Gs(ojk)2  T2 [Q1(u>k,  t,  O1)  Q2(wk,  t,  £2)  -  Q1>2(wk,  T,  D1,  D2)  • 
Qlj2((JJk»  T,  D1,  £2)]  +  Gs(wk)  x[G  (a^)  Qj »  T>  21) 

*  G^f^)  Q2<V  t,  I2)!  *  G^)  en2(a,k)  .  (3-32A) 

We  can  now  express  the  conditional  density  function  P(Y/D  ,  £  ) 
(Equation  3-10)  by  substituting  Equations  3-32  and  3-32A  into 
Equation  3-10  via  Equations  3-10A  and  3-10B. 

N 

Pd/O1,  O2)  =  J~J  [2tt|  Z(uok)  i  j"1  exp{  t-C^T(cok)  Z'1^)  CfUk)]1/2} 

k=-N 

K  n  (3-33) 


The  standard  procedure  is  to  take  the  natural  log  of  the  condi¬ 
tional  density  function  (Equation  3-33).  We  shall  also  let  t  effec¬ 
tively  approach  infinity  while  maintaining  _j.  =  2-rrk/ x  as  a  constant. 


Note  the  following  limits: 


“k  ! 
wN  = 

"-N 

Lim 

T-v® 

Lim 

t-h» 

where  X^i 
sensor  1 
sensor  2 

The 


■  2iU  u. 

T 

(3-34A) 

.  M-  B 

T 

(3-34B) 

.  -2irN  , 

=  -  +  | 

T 

3 

(3-34C) 

cf  ^  (“k )  = 

(3-34D) 

T  cf2K> 

((i^)  “  ^2 (co) 

(3-34E) 

u)  is  the 

Fourier  transform  of  received  signal 

x^(t)  at 

and  XoU) 

is  the  Fourier  transform  of  received 

signal  x2(t)  at 

~1 

ML  estimator  are  the  parameters  0  ,  £  which  maximize  the 


following  expression  resulting  from  taking  the  natural  log  of 
Equation  3-32A  and  converting  summation  to  their  appropriate  integrals. 


3-16 


§  f  ln(^)  dco 

J.R 

(1) 


-  |  f  1  n I  z((jj)  i  du> 

-B 


(2) 


D 

_1_  f 

4tt  / 
l  -B 


[ | X^Cuj)  1  Z2  2^)  +  ix2^a)^  2l,l(“)]dw 

*  FzT^TI 


(3) 


D 

h  f  xi(a))x2 

-B 


(«) 


z1.2(ai)  +  z2.1(a)),J 
-r\zm — |dw 


(4) 


(3-35) 


where  z-  •  are  the  elements  of  matrix  z(u)  given  in  Equation  3-31. 


Note  expression  3-35  is  symbolically  the  same  as  that  arrived  at  by 
2 

Knapp  and  Carter. 


Let  us  now  examine  each  of  the  four  terms  in  expression  3-35. 


~1  ~2 

The  first  term  is  independent  of  D  arid  £  and  therefore  can  be 


ignored. 


The  second  term  is  a  function  of  the  determinant  |Z(co)|.  This 
determinant  as  can  be  seen  frcm  Equation  3-32A  is  a  function  of  0  and 
D2  but  is  also  a  function  of  the  noise  power  spectrums  G  (oi)  anc 


•..s  ‘.‘.V.Vl 


-  '  . 


UTf, 


G  (to).  Although  it  is  not  clear  at  this  point,  for  reasonable  target 
n2 

c  ~i  -2 

motion,  the  primary  effect  of  non-zero  £  and  2  is  a  small  distortion 

in  magnitude  of  the  received  auto  and  cross-power  spectrums.  Certainly 

~2 

at  low  SNR,  variations  in  term  2  due  to  £  and  will  be  unobservable 

even  if  a  significant  dependence  exists.  Note  that  for  no  target  motion 
-1  -2 

[D  =0  =0]  and  our  flat  signal  and  noise  power  spectrum  assumption 
that  | z (oj) |  is  a  constant  for  all  frequencies  of  interest.  For  now,  we 
shall  ignore  term  2  and  examine  term  3. 

~1  ~2 

Term  3  is  dependent  on  D  and  0  but  is  also  dependent  on  the  noise 

power  spectrum  Gn  (to)  and  (co).  Examining  term  3,  it  symbolically 
nl  n2 

appears  to  be  the  sum  of  cross-matching  the  measured  auto  power  spec¬ 
trums  for  sensors  1  and  2  against  what  can  be  interpreted  as  the 
expected  auto  power  spectrums  for  sensors  1  and  2.  For  small  expected 

auto  power  spectrum  magnitude  variations,  term  3  will  be  weakly  depend- 
-1  -2 

ent  on  D  and  D  and  certainly  unobservable  at  low  SNR.  Therefore,  we 
shall  ignore  term  3. 

~1  ~2 

Term  4  is  dependent  on  D.  and  D  and  is  not  significantly  dependent 

on  the  noise  power  spectral  densities  G_  (a)  and  G_  (to).  In  fact,  as  we 

nl  n2 

shall  see  later,  term  4  is  primarily  dependent  on  the  difference  param¬ 
eter  AD,  which  is  exactly  what  we  want.  Therefore,  the  estimate  5^  and 
-2 

jD  which  maximizes  term  4  shall  be  our  primary  interest.  Knapp  and 
2 

Carter  reached  a  similar  conclusion  for  a  first  approximation. 


As  noted  earlier,  the  determinant  jz(u)j  is  weakly  dependent  on 
-2 

and  £  and  this  dependency  is  certainly  unobservable  at  low  SNR.  In 
addition,  our  flat  noise  and  power  spectrum  assumptions  makes  |Z(oj)| 
effectively  constant  as  a  function  of  co.  For  this  reason,  we  shall 
write  the  ML  estimator  for  term  4  without  |z(gj)|. 


T-z  k  / 

d\[T  *'-b 


X ^ ( cj )  X2(u)) 


2^  2(to) 


zl,2(oj) 


dcu 


(3-36) 


Using  Parsevals  theorem,  we  can  view  Equation  3-36  as  equivalent 

-1  - 2 

to  maximizing  with  respect  to  D  and  D  the  integral  of  the  product  of 
the  measured  cross-correlation  function  with  the  expected  cross¬ 
correlation  function. 


MAX 


/ 

-CD 


R(t)  R(t)  dt 


(3-37) 


where  R(t)  is  the  measured  cross-correlation  function  between  sensors  1 
and  2;  and  R(t),  the  expected  correlation  function,  is  given  by  the 
following  inverse  Fourier  transform. 


s1>2(u>)  ♦  z1>2( 


dcu 


(3-38) 


Equation  3-38  shall  now  be  referred  to  as  the  mismatch  function. 


We  now  have  a  general  form  for  a  practical  approximate  ML  estimator 
for  the  unknown  target  time  delay  difference  trajectory.  In  the  next 


|OI 


section,  we  shall  examine  a  more  pliable  target  motion  model  and  obtain 
a  precise  expression  for  the  mismatch  function. 

3.4  THE  MAXIMUM  LIKELIHOOD  ESTIMATOR  FOR  PROBLEM  2 


From  a  dynamic  and  a  practical  point  of  view,  it  is  not  desirable 

to  wait  an  entire  x  seconds  to  complete  our  measurements  and  make 

decisions.  Let  us  assume  that  we  actually  make  M  measurements  of  T 

seconds  each  over  the  x  second  interval  [i.e.,  M  T  =  x] .  We  shall  also 

assume  that  the  unknown  time  delay  difference  trajectory  is  a  piecewise 

linear  function  as  described  in  Equation  2-6,  with  a  linear  segment  for 

1  2 

each  of  the  M  time  segments.  The  unknown  parameter  vectors  D  and  D 

1  2 

now  represent  the  time  delay  components  dfe,  d^  at  the  end  points  of  each 

of  the  k  =  1,  2,  3,  ...  M  piecewise  linear  segments.  In  addition,  the 

12  12 
initial  time  delays  dg,  dg  may  be  considered  part  of  D  and  D  .  As  an 

example,  vector  D*  components  are  shown  below. 


M1  H1 

[dg,  d^ 


(3-39) 


We  will  now  obtain  an  approximation  of  Equation  3-37  with  our  new 
time  delay  trajectory  model.  Since  the  cross  correlation  is  a  linear 
process,  the  measured  cross-correlation  function  R(t)  for  the  entire 
interval  x  can  be  viewed  as  the  linear  sum  of  the  M  subintervals 


measurements  Ar. 


[•] 


For  the  same  reason  the  expected  cross-correlation  function  (mis¬ 
match  function)  can  be  viewed  as  the  sum  of  each  of  the  M  subinterval 
expected  cross-correlation  functions  Ar  .j . 

M 

.  (3-4D 

i=l 


Substituting  Equations  3-40  and  3-41  into  3-37  yields  the  follow¬ 
ing  result: 


MAX 
1  2 

D  ,D  1-1 


M  ® 

s/ 


Ar i ( t )  Ar.(t)  dt 


M 

i=l  j=i+l 


M  00 

£./ 


Ar i  ( t )  Ar j(t)dt 


(3-42) 


Time  segments  Ar,  and  Ar  •  for  i  f  j  are  separated  by  at  least  T 

1  J 

seconds.  We  shall  assume  that  the  contribution  to  the  overall  correla 
tion  function  due  to  these  terms  is  minimal  and  thus  Equation  3-42  can 
be  approximated  by: 


MAX 

D1^2 


1*1  ® 

Y J  4ri (t)  Ar~(t)dt 
i=l  -<» 


(3-43) 


Note  that  the  contribution  of  /  Ar._^(t)  Ar^tjdt  will  certainly 


be  greater  than  J"  Ar.  ^(t;)  Ar^(t)dt  which  will  be  greater  than 

•cd 

00 

f  Ar^U)  Ar.(t)dt  for  most  actual  target  trajectories .  This  assump- 

-oo 

tion  will  be  of  use  later  when  determining  a  suitable  a  priori  informa¬ 


tion  penalty  function. 


Equation  3-43  is  attractive.  The  integral  in  Equation  3-43  can  be 
replaced  by  a  summation  yielding  a  form  conducive  for  batch  techniques. 
In  addition,  the  expected  cross  correlation  or  mismatch  function  Ar^t) 
is  now  defined  for  a  single  piecewise  linear  segment.  We  shall  take 
advantage  of  this  development  to  derive  the  detail  form  of  the  mismatch 
function. 


Let  us  now  rewrite  our  receive  signal  model  to  reflect  the  above 
assumptions.  For  the  first  segment  we  have  the  following: 


X]_(t)  =  s ( t  +  ^(D1,  t))  +  n1(t ) 


x2(t)  =  s(t  +  <p2(D  ,  t))  +  03(0 


where 


and 


MAX  |  As  1 1  and  MAXjASj  <<  1 


(3-44) 


m 


t)  =  dj  +  Asj  t,  0  <  t  <  T 

(3-44A) 

1  1 

t)  »  d*  +  As*  t,  0  <  t  <  T 

( 3-44B) 

v  3-44C ) 
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3-22 

Although  not  required,  we  shall  assume  for  convenience  that  the 

1  2 

initial  received  time  delays  dg  and  dg  are  known  for  the  first  subin¬ 
terval.  For  all  subsequent  subintervals,  the  initial  time  delays 
d^  and  d^  will  be  determined  from  the  previous  subinterval  (i.e., 
d^  =  d^  1  +  As^  ^  T  and  d^  =  df_1  +  As^  T).  Since  all  the  piecewise 
linear  segment  are  identical,  we  shall  proceed  with  our  analysis  using 
the  first  segment  as  a  model. 

The  first  step  in  calculating  the  precise  form  of  the  mismatch 
function  is  to  calculate  the  precise  form  of  the  non-stationary  Fourier 
coefficients  a^l)  and  a^(2)  for  a  linear  moving  target. 


’  .  - 


Let  us  define: 

82  =  1  +  Asj  =  1 
B2  =  1  +  As^  =  1 
and  note  that 

-  1  T  -  dJ 
4>!(D  ,  y)  =  ~gp 

-  ,.2  ,  T  -  dg 

*  y'}  =  r2 

T  =  t  . 


(3-45A) 

(3-45B) 
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(3-45C) 

(3-45D) 


(3-45E) 
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We  can  now  substitute  Equations  3-45A  through  3-45E  into  the 
expression  for  1^  (Equation  3-19A)  yielding  the  following  result. 


h  - 


■u 


J  -a.A ( Y-d J ) /Bx ] 


dy 


(3-46) 


Evaluating  the  integral  in  Equation  3-46  and  with  the  appropriate 
algebraic  manipulation  yields  the  following  expression  for  I,. 


I 
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„JtVo 


j  ( 3-i o)_ -u),  )  (T/2) 

Sinc[(S1  -  cok)  (T/2)  ]  e  1  1  K 


(3-47) 


Note  Sinc(x)  is  defined  as  sin(x)/x. 


Similar  manipulation  can  be  performed  yielding  an  identical  form 
for  Ig. 


>zdo 


Sine  [(3^  -  .juk )  ( T/2 )  ]  e 


J(32a)l':j!<)(T/^) 


(3-43) 


Equations  3-19  and  3-20  are  still  valid  for  the  non-stationary 
Fourier  coefficients  o^(l)  and  a.,  (2),  respecti ve ly ,  but  now  use  the 
above  expressions  for  1^  and  I?  which  are  based  on  the  linear  time  delay 
trajectory  model . 


The  next  step  in  our  calculation  is  to  derive  a  precise  expression 

*T 

for  the  covariance  matrix  E[j3f(wk),  (tuk ) ]  (Equation  3-12).  -The  (1,1) 
element  in  this  matrix  is  given  by  Equation  3-25.  If  we  substitute 
Equation  3-47  into  Equation  3-25  and  noting  that  T  =  t  we  have  the 
fol lowing. 

B 

E[cf  (u)k),  cf  (ujk)]  =  -A  f  6s(w)  Lim  (Sine2  [(S^oj  -  wk)^]  )daj 

1  1  "  g  T-^30 

+  Gnx(a)k)/T  ‘  (3-45) 

~  G1  ^i 

Now  let  uj  =  —  to  and  dw  =  j-  dej  and  substitute  into  Equation  3-49 


E(Cf1(uk)’  =  2%  f 


BB/2 


-86/  2 


G A=^)  Lim(Sinc4 [(w 
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(3-50) 


As  we  let  T  approach  infinity,  the  limit  of  the  Sine  (■)  approaches 


a  delta  function4  as  shown  below. 


<*>L 


2tt 


Lim(Sinc  [ (oj  -  -^)Tj)  -  DeltaU  -  -4 


(3-51) 


If  we  now  assume  that  G$(<jj)  is  relatively  constant  over  the  inter 

QV  Q  (jjl  O 

val  — y»  +  >  we  can  substitute  Equation  3-51  into  Equa¬ 

tion  3-50  and  evaluate  the  integral  to  yield  the  following. 


E [cf x (^k ) ’  cf1^a)k^  =  b/t  Gs^  +  T  <  wk  <  B 

(3-52) 

The  above  result  is  pleasing.  It  states  that  the  doppler  shift  due  to 

linear  target  motion  results  in  either  an  expansion  or  contraction  of 

the  stationary  target  auto  power  spectrum  with  an  appropriate  scale 

factor  to  maintain  constant  total  power.  This  result  is  in  agreement 

2 

with  Knapp  and  Carter. 

Certainly  for  reasonable  signal  power  spectrums  passed  through  a 
prefilter,  the  slight  level  change' (1/8^)  in  the  received  signal  power 
spectrum  will  be  completely  unobservable  with  respect  to  other  natural 
variations  in  received  signal  power  levels  or  noise  power  levels.  Thi 
result  supports  our  decision  to  ignore  term  2  in  expression  3-35. 

An  equation  similar  to  Equation  3-52  can  be  derived  for  the  (2,2) 

★T 

element  of  the  covariance  matrix  E(Cf(cj|<),  C^  (u^)]  which  is  given 
below. 
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The  critical  (1,2)  element  of  covariance  matrix  E[Cf(wk),  Cf  (wk)] 
is  more  difficult  to  analyze.  Substituting  Equations  3-47  and  3-48 
into  Equation  3-28  yields  the  following  key  integral. 


D 

E[cf  (o^j,  cf  (ojk)J  =  ^  f  Gs(cu)  Lim(Sinc  [(e^oj  -  cuk)|-]  * 
1  2  J  a  t** 


T  jaj(Ad0+AB(T/2)) 

Sine  [(82  w  "  tu^^le  )du>  . 


Note  Adg  =  dg  -  dg  and  A8  =  8^  -  S2 


(3-54) 


NgiJ  developed  an  excellent  direct  analytical  approximation  to  the 
integral  of  Equation  3-54  which  is  shown  in  slightly  modified  form  in 
Appendix  A.  This  direct  analytical  result  shall  be  shown  after  the 
following  intuitive  approach  which  was  initially  used  to  approximate 
the  integral  of  Equation  3-54. 


The  key  concept  comes  from  examining  the  method  actually  used  in 
practice  to  calculate  a  correlation  function  (correlogram)  over  a 
period  of  time  T.  The  actual  technique  is  to  average  consecutive  sums 
of  m  smaller  correlograms  each  of  duration  At  seconds.  (In  practice,  At 
is  determined  to  obtain  a  sufficient  bandwidth  time  product;  however, 
for  reasonable  bandwidth.  At  can  be  quite  small.)  Therefore  we  have: 


yv 


The  integrand  of  Equation  3-54  represents  a  process  observed  over 
T  seconds.  Let  us  redefine  the  integrand  of  Equation  3-54  as  the  aver¬ 
age  of  the  sum  of  m  processes  observed  consecutively  every  At  seconds. 

We  shall  assume  that  for  each  process  the  time  delay  between  sensors  is 
constant  for  the  given  At  second  interval.  This  assumption  implies  that 
for  each  subinterval : 

=  82  =  for  convenience  MAXfg^,  g^j  =  g  (3-56) 

Note:  At  this  point  there  is  no  firm  guidance  on  what  to  choose 
for  the  constant  value  of  B  although  it  is  obvious  that  ‘ 

Min(6^,  S2)  <  3  <  MAX(g^,  g^).  We  choose  Equation  3-56  to  be  compat¬ 
ible  with  the  result  obtained  by  Ng.^  For  our  assumptions  on  g^  and  32 
the  above  is  a  mute  point. 

However,  the  constant  delay  for  each  At  subinterval  is  given  by: 

del  ay i  =  AdQ  +  i  At  A8,  i  =  1,  2,  3  ...  m  .  (3-57) 

This  model  represents  the  linear  time  delay  difference  trajectory 
by  a  staircase  approximation.  The  validity  of  this  model  is  dependent 
on  8^  =  S>2  and  sufficient  bandwidth. 


Note  that  the  Sinct(*)  functions  in  Equation  3-54  now  have  ident 


i cal  arguments  which  will  remain  unchanged  for  all  m  subintervals. 

We  shall  now  let  i  =  3  cd/2  and  dcu  =  3/2  dm.  Substituting  into 
Equation  3-59  yields: 
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Sinc^[(a)  -  -f)At]  ^ 
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(2w/6)(AdQ+iAtA3) 
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(3-60) 


Now  if  we  let  At  become  relatively  large,  the  exponential  terms  are 

p 

not  unduly  affected  but  the  Sine  (•)  term  approaches  the  following 
limit. 


Lim  Sinc^[(u  -  -|)At]  -*  ^  Delta(al  -  — =•)  .  (3-61) 

At-*®  ~ 


For  the  exponential  terms,  if  we  now  let  m  become  large,  we  have  by 
the  definition  of  an  integral  and  Equation  3-55  the  following  limit: 


v/v: 


Lim 

m+a> 


j(2w/8)(Ad0+(i/m)A8T) 


j ( 2cj/ 8 )  (Adn+xA8T) 
e  u  dx 


c.  ru)A8T  j(2w/8)  (Ad0+(A8/2)T) 
■  S  i  nc  ]  e 


13-62) 


Substituting  Equations  3-61  and  3-62  into  Equation  3-60  we  have  the 
following  result. 
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If  ?  n  ~AftT  j(2aj/8)(AdQ+(A8/2)T) 

8?  /  Gs(  8  ]  De1ta(oj  '  ~7]  Sinc[^  e 
-SB 


doj 
(3-63) 


If  we  now  assume  that  Gs(uj)  is  relatively  constant  over  the  inter- 


val  (-p  -  2~  +  ^ ] ,  Equation  3-63  evaluates  to: 


Eicf1K>- K>i  =irGs(r)  S1"c( 


coj^  A6  T 
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j(o)k/3)(Ad0+(A8/2)T) 


-SB  <  uik  <_  8B 


(3-64) 


The  direct  evaluation  of  Equation  3-54  by  Ng^  yields  the  identical 


result  with  similar  assumptions.'  (See  Appendix  A.) 


Equation  3-64  yields  intuitively  pleasing  results.  Similar  to  the 
case  for  tne  received  auto  power  spectrum.  Equation  3-64  predicts  a 
frequency  expansion  or  contraction  with  an  appropriate  scale  factor  to 
maintain  constant  total  power.  This  effect  is  primarily  due  to  doppler 
and  is  inversely  proportional  to  S  =  Max(3]_.  Sj).  Equation  3-64  also 
predicts  a  loss  of  coherency  primarily  due  to  tangential  target  motion 
represented  by  AS-  Finally,  Equation  3-64  predicts  an  effective  time 
delay  shift  to  the  midpoint  of  the  observation  interval  represented  by 
Ad0  +  T*  The  Sinc(-)  factor  in  Equation  3-64  is  the  significant 
difference  between  our  result  and  the  result  obtained  by  Knapp  and 
Carter^ 

Similar  to  past  analysis,  we  have  that  3=1  and  thus  Equa¬ 
tion  3-64  can  be  rewritten  to  a  much  more  convenient  form. 


E[cf  (4^),  cf^(ujk)3 


.((A.)  juv(Adn+(AS/2)T) 

Sinc[<^  AST/2]  e  K  u 


Assuming  Adg  is  known,  Equation  3-65  is  only  dependent  on  one  param¬ 
eter,  AS,  or  equivalently,  AST. 

Equation  3-65  via  Equation  3-31  can  be  substituted  into  Equa¬ 


tion  3-38  to  yield  an  expression  for  the  mismatch  function. 


R(t)  =  ~  I  Gs(u)  Sinc(o)  T)  cos[u(t  +  [AdQ  +  f-  T])]dw 

J  D 


(3-66) 


We  shall  now  make  use  of  the  flat  band-limited  assumptions  on  Gs(ui) 
given  in  Equation  2-3  and  repeated  below  for  convenience. 


Gs(f)  ■ 


0  ,  elsewhere. 


(3-67) 


Let  =  2itf^  and  ^  =  2 7rf 2 


OJj  +  (Jig 

Bw  =  (^2  -  and  u>c  =  - ^ — “ 


and  for  now  assume  that  the  total  received  signal  power  A  is  1. 


Substituting  the  above  assumptions  on  Gs(u)  into  Equation  3-66 
yields  the  following  after  some  algebra. 


R(t)  -  h 


i-  r 

Bw  J 

'll  - 


ujc+8w/2 


Sinc(co  M  j)  cos (ui[t  +  (AdQ  +  f-  T) ]  )d 


ajr-8w/2 


This  integral  is  analytical ly  approximated  in  Appendix  3,  yielding 
an  analytical  equation  for  the  mismatch  function  which  depends  on  the 
following: 

1.  t:  the  correlation  lag  point. 

2.  Adg:  The  initial  time  delay  difference  between  the  two 

sensors. 

3.  ABT:  The  total  change  in  the  time  delay  difference 

between  the  two  sensors  during  the  observation  time 
interval . 

Since  Adg,  T  are  assumed  known,  then  the  only  quantity  to  maximize  Equa 
tion  3-37  is  A6  or,  more  conveniently,  A0T. 

At  this  point  we  shall  summarize  what  we  have  derived.  If  we 
assume  that  the  unknown  signal  time  delay  difference  trajectory  is 
linear  over  an  observation  interval  of  T  seconds  with  a  known  initial 
time  delay  difference,  we  can  determine  an  approximate  ML  estimate  for 
the  unknown  slope  of  the  time  delay  difference  trajectory  by  maximizing 
the  product  between  the  measured  cross  correlation  function  and  the 
expected  cross-correlation  (mismatch)  function  with  respect  to  the 
unknown  slope.  Also,  we  have  determined  that  an  approximate  ML  esti¬ 
mator  for  an  unknown  piecewise  linear  time  delay  difference  trajectory 
(Problem  2)  is  to  maximize  the  sum  of  the  individual  but  continuous 
linear  segment  correlations  with  respect  to  their  unknown  slopes.  In 
addition,  we  have  derived  an  approximate  analytical  equation  for  the 
expected  cross-correlation  function  (mismatch  function)  assuming  a 


linear  time  delay  difference  trajectory,  relatively  slow  target  motion, 
and  sufficient  signal  bandwidth. 


3.5  LEAST  MEAN  SQUARE  FIT  APPROXIMATION 


At  this  point  it  is  productive  to  reformulate  our  estimator  in  a 
similar  but  more  practical  form.  Our  input  measurements  over  an  inter¬ 
val  of  time  T  will  take  on  the  form  of  normalized  noisy  discrete  cross¬ 
correlation  lag  points  designated  by  the  vector  R,  where  the  components 
r.  of  the  vector  R  are  measurements  of  the  cross  correlation  at  discrete 
times  tj  -  i  tau,  i  =  -N,  -N  +  1,  . . . ,  0,  . . .  N  and  tau  is  a  constant 
time  delay.  The  normalization  factor  shall  be  an  estimate  of  the  total 
received  power  at  either  sensor  1  or  2  (signal  power  (A$)  +  noise  power 
(Ns)).  Therefore  at  tQ  =  0, 

rQ  =  As/(As  +  Ns)  (3-69) 

and  for  stationary  conditions  (no  target  motion)  we  have 


Lim  rQ  =  A$/(As  +  N$)  .  (3-69A) 


In  addition,  we  shall  at  least  initially  assume  that  the  noise 
component  of  the  lag  measurements  r^'s  is  correlated  and  has  a  known 
constant  correlation  matrix  W.  This  assumption  corresponds  to  filter¬ 
ing  and  over-sampling  our  original  noise  processes  r^U)  and  r^U) 
(Chapter  2).  We  shall  show  that  the  correlated  noise  problem  can  be 
transformed  into  an  independent  noise  problem. 


V-V-V' 


»■- 


m  JL — uuk. 


J.  JLW.U  ’ 


r-v-'V  V-/ 


Our  new  estimator,  shown  below,  will  minimize  the  squared  error 


between  the  measured  cross-correlation  lag  data  R  and  a  normalized 
expected  cross-correlation  lag  data  A  *  £  with  respect  to  the  unknown 
time  delay  difference  slope  parameter  A8T  and  the  normalization 
factor  An. 

MIN  [R  -  A  C]T  w"1  [R  -  A.  C]  (3-70) 

A6T,A 

where  the  components  c  •  of  the  vector  £  are  the  mismatch  Equation  3-68 
evaluated  at  times  t..  =  i  tau,  i  =  -N,  -N  +  1,  . . . ,  0,  . . .  N. 


Note  that  the  unknown  normalization  factor  An  is  an  estimate  of 
the  signal  power  divided  by  the  total  receive  power. 


An  = 


A  +  N 


(3-71) 


If  w  is  the  identity  matrix,  our  measured  data  is  uncorrelated  and 

we  have  a  straightforward  least  squares  problem.  If  W  is,  as  we  assume 

a  constant  symmetric  positive  definite  matrix,  we  can  transform 

Equation  3-70  into  a  straightf orward  least  squares  problem.  This 
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technique  is  called  pre-whitening  and  is  described  by  Whalen  and 
briefly  outlined  below. 


3-36 


If  W  is  positive  definite,  we  can  always^  factor  it  into  lower 
(L)  and  upper  (u)  triangular  matrices  which  are  transposes  of  each 
other  (L  =  U).  These  triangular  matrices  can  be  easily  inverted 
to  yield  the  following  equation 

W"1  =  L-1T  L-1  .  (3-72) 

If  Equation  3-72  is  substituted  into  Equation  3-70  along  with  the 
following  substitutions, 

R  =  L'1  R 

C  =  L'1  C  , 

we  arrive  at  the  following  simple  least  squares  problem. 

MIN  [R  -  An  C]T  [R  -  A  C]  .  (3-75) 

A6T,An  n  -  n  - 


(3-73) 

(3-74) 


*  The  factorization  process  is  a  form  of  a  Choleski  decomposition. 
In  practice,  a  small  positive  factor  must  be  added  to  the  main 
diagonal  terms  in  order  to  ensure  numerical  positive  definiteness 
and  proper  factorization. 


'  .  V  v  •  .  <.  **. 
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3-37 


The  effect  of  Equation  3-73  is  to  uncorrelate  or  whiten  the  noise 
components  of  the  measurement  vector  Rl.  Since  the  whitening  process  of 
Equation  3-73  affects  the  underlying  measured  signal  cross-correlation, 
Equation  3-74  corrects  the  expected  signal  cross-correlation  model. 


From  this  point  onwards,  we  shall  talk  to  Equation  3-75  and  its 
assumed  assumption  of  uncorrelated  noise.  If  correlated  noise  is 
present,  we  have  the  option  of  implementing  Equations  3-73  and  3-74  to 
transform  the  problem  to  the  case  of  uncorrelated  noise. 


To  minimize  Equation  3-75  with  respect  to  An,  we  shall  take  the 
partial  derivative  of  Equation  3-69  with  respect  to  An  and  set  the 
result  equal  to  0.  The  resulting  estimate  for  An  is  given  below. 


(3-76) 


Substituting  Equation  3-76  into  Equation  3-75  yields  after 
expansion: 


MIN  RT  R 
AST 


(CT  R)2 
CT  C 


(3-77) 


Only  the  second  term  in  Equation  3-77  depends  on  the  parameter 
AST.  Unfortunately,  the  solution  is  not  equivalent  to  maximizing  the 
lag  product  of  the  matching  function  C  with  the  measured  cross¬ 
correlation  data  R.  This  is  true  because  the  dot  product  C  is  not 
constant  for  all  choices  of  AST. 


However,  we  can  modify  the  definition  of  our  matching  function  to  the 


fol lowing . 


Matching  function  = 


(3-78) 


where 


A  =  A  \C  C  . 
n  n  " —  — 


Equation  3-75  can  be  rewritten  to  the  following  form. 


(3-7SA) 


MIN  [R  -  A  C/  V?  C] T  [R  -  A  C/  Jc  C] 
A0T,A  n  1  •  n  * 


(3-79) 


If  we  take  the  derivative  of  Equation  3-79  and  set  the  result  equal 
to  zero,  we  obtain  the  following  expression  for  the  estimate  of  A  . 


(3-80) 


Substituting  Equation  3-80  into  Equation  3-79  and  expanding  yields 
the  fol  lowing  result. 


T  (CT  R)2 

MIN  R  R - j - 

A0T  C 1  C 


( 3-si ; 


This  is  the  same  expression  as  Equation  3-77,  but  now  the  second 
term  can  be  rewritten  as  the  following. 


Maximizing  Equation  3-82  is  equivalent*  to  maximizing  the  follow¬ 


ing  expression. 


Equation  3-83  is  exactly  what  we  want;  maximizing  the  lag  product 
between  normalized  matching  function  and  the  normalized  measured 
cross-correlation. 

This  result  can  easily  be  extended  to  the  case  of  a  piecewise 
linear  time  delay  difference  trajectory  over  a  block  of  M  data  segments 
With  Equation  3-43  in  mind,  define  the  measured  cross-correlation 
vector  for  time  segment  j  as  R,.  AST  now  uecomes  a  vector  with  a 

J 

component  for  each  unknown  piecewise  linear  slope.  The  approximate  ML 
estimator  can  now  be  written  as: 

M  _ _ 

MAX  R'.  C/  VCT  C  .  (3-84) 

Note  that  C  has  an  implied  dependency  on  the  time  segment  j  since 
eacn  time  segment  has  a  unique  slope  parameter  A3-;  i . e . ,  C(Afi.T). 

J  J 


*  Note  the  above  expression  disregards  negative  maxi  mums,  but  these 
are  not  of  interest. 


■0 


Equation  3-84  can  be  rewritten  in  an  equivalent  form. 


(3-85) 


where  A  is  the  estimated  power  normalization  factor  for  the  jth  seg- 

nj 

ment.  With  Equation  3-71  in  mind,  the  SNR  estimate  for  the  jth  interval 


is  given  by: 


(3-86) 


The  average  SNR  estimate  for  the  entire  block  of  data  is  given  by: 


M 

SNRAVG  *  S  2 


j=l 


(3-87) 


In  summary.  Equations  3-84,  3-85,  and  3-87  are  the  key  equations 
for  the  practical  ML  estimation  portion  of  the  MAP  estimator 
(Equation  3-3).  Equations  3-84,  3-85,  and  3-87  assume  a  piecewise 
linear  target  time  delay  difference  trajectory. 


3.6  SPECTRAL  PARAMETER  ESTIMATION 


In  the  previous  sections  we  have  assumed  that  the  signal  power 
spectrum  G$(w)  can  be  modeled  by  a  flat  bandpass  spectrum  given  by 
Equation  2-3  and  repeated  below  for  convenience. 
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(3-88) 


Note  oj  =  2rf. 

It  would  not  be  a  serious  violation  of  previous  work  if  we  let  the 
signal  power  spectrum  Gs(cj)  have  one  more  degree  of  freedom  by  allowing 
an  unknown  spectral  slope  parameter  y. 


Let  us  model  the  signal  power  spectrum  G  (w)  by  a  linear  model 
shown  in  Figure  3-1  for  a  one-sided  spectrum. 


Figure  3-1.  Linear  Model  for  the  Signal  Power  Spectrum 


Let  us  define  f 


(3-90) 


Then  the  one-sided  power  spectrum  G  (•)  can  be  defined  for  any 
frequency  as  the  following. 


.  - . *■. •-«  v. r. 


of  Equation  3-68  becomes  the  following. 
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cos(a)[t.  +  Adn  +  ^p])  sin(a>^)du 
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ci).+Bw/2 
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cos(a)[t.  +  Ad q  +^])  sinc(co  ^jp)do 


(3-94) 


The  first  and  third  integrals  are  effectively  equivalent  to 
Equation  3-68.  All  integrals  are  analytically  or  approximately  ana¬ 
lytically  evaluated  in  Appendix  B.  Equation  3-94  yields  a  new  form  for 
the  mismatch  function  represented  by  the  vector  C.  If  we  designate  the 
solution  of  the  first  integral  at  lag  time  ti  as  u.  and  the  solution  of 
the  second  and  third  integrals  at  lag  time  t.  as  v ^  we  can  write  the 
following  expression  for  the  new  components  c^  of  the  vector  C. 


ci  =  ui  +  ^vi 


(3-95) 


C  =  U  +  y  V  . 


(3-95A) 


Now  we  substitute  Equation  3-95A  into  Equation  3-33  to  obtain  a 


revised  form  of  the  ML  estimator. 


MAX  (U  +  yV)T  RT/((U  +  yV)T  (U  +  yV))1/2  . 

A3T 


(3-96) 


Equation  3-98  can  be  maximized  numerically  for  A6T,  yielding  an 
estimate  for  AST  and  also  indirectly  via  Equation  3-97  an  estimate  for 
y.  However,  Equation  3-98  is  significantly  more  complicated  than 
Equation  3-83  and  contains  implicitly  an  estimate  of  the  rather 
unstable  slope  parameter  y. 

A  second  more  practical  technique  is  the  following. 

Assume  that  the  slope  is  zero;  i.e.,  y  =  0.  This  assumption  via 
Equation  3-96  yields  our  original  ML  estimator  given  in  Equation  3-83. 
Note  that  under  the  single  target  assumption  the  maximum  of  Equa¬ 
tion  3-83  will  occur,  on  average,  when  the  expected  peak  of  the  measured 
cross-correlation  R  lines  up  with  the  symmetrical  peak  of  the  expected 
cross-correlation  vector  C.  The  effect  of  the  spectral  slope  parameter 
(y)  is  to  broaden  or  narrow  the  cross-correlation  function.  It  does  not 
affect  the  location  of  the  peak.  Therefore,  Equation  3-96  can  be  optim¬ 
ized  independently  of  the  slope  parameter  (y).  Once  this  is  completed. 
Equation  3-97  can  be  used  to  obtain  an  estimate  of  the  spectral  slope 
parameter  y. 

3.7  DERIVATION  OF  THE  A  PRIORI  KNOWLEDGE  PENALTY  FUNCTION 

We  have  obtained  an  expression  for  the  ML  portion  of  the  MAP 
Equation  3-3.  We  shall  now  obtain  an  expression  for  the  3  priori 
knowledge  portion  of  Equation  3-3. 


In  essence,  the  a  priori  knowledge  is  what  ties  Problem  2  and 
Problem  1  together.  Our  ML  estimator  assumes  a  piecewise  linear  time 
delay  difference  trajectory  and  yields  an  estimate  of  the  time  delay 
difference  for  each  piecewise  linear  segment  (Problem  2).  However,  we 
assume  a  priori  in  Problem  1  (Equations  3-4,  3-5)  that  the  actual  time 
delay  difference  trajectory  is  quadratic;  i.e., 

Adi  =  Ad2  t^  +  Adj^  ti  +  Adg  .  (3-99) 

Note  that  estimates  of  Ad ^  from  the  ML  estimator  (Equation  3-84) 
are  given  by 

i 

Adi  =  AdQ  +  ^  ABjT  (3-100) 

j=l 

where  AdQ  is  the  estimated  or  given  time  delay  difference  at  the  begin¬ 
ning  of  the  trajectory  and  A8T  are  the  estimates  obtained  from 
Equation  3-84. 

We  shall  assume  that  the  difference  between  the  ML  estimates  Ad^ 

and  the  truth  Ad^  are  independent  Gaussian  random  variables  with  equal 
2 

variances*  a  .  Therefore,  the  a  priori  density  function  P^q(AD)  is 
given  by  the  following  expression. 

*  We  have  assumed  that  all  the  components  of  AD  are  equal  variance  and 
are  uncorrelated.  Neither  assumption  is  entirely  true,  so  theoret¬ 
ically  there  is  a  covariance  matrix  Q  relating  the  components  of 

ADi.  Equations  3-102,  3-103,  and  3-104  should  therefore  be  modified 

by  the  components  of  the  Q  matrix.  However,  the  ML  estimates  AD  and 
the  a  priori  information  interact  in  a  complex  manner  whic.i  obscures 
this  relationship.  Simulations  have  shown  small  gain,  if  any,  for 
weighting  matrixes  Q  other  than  the  Identity. 


^■ry= 


“T  i  -[Ad^Ad^+Ad^+Adp)]2 


(3-101) 


1-1  '2w 


PAI)(^  •  M  ln<-=>  -  52  Z  C^i  •  (AdV?  +  *Vi 

—  /2to  i-1 


+  Ad0)]^  . 

(3-102) 


Only  term  2  of  Equation  3-102  is  dependent  on  AD  or  equivalently 
AgT.  Therefore,  finally,  the  MAP  Equation  3-3  can  be  written  as  a 
combination  of  Equations  3-102  and  3-84. 


M  C 


MAX  -J  -  K  [Ad.  -  (Ad,t2  +  Ad,  t .  + 

arT  L-J  J  2j  1  j 

—  j*i  Vcr  C  j-1 


Ad  n )  j  (3-103) 


We  can  rewrite  the  second  summation  in  Equation  3-103  in  the 
following  standard  matrix  format. 


yt  [Ad j  -  (Ad 2  t2  +  Ad1  t-  +  Ad0)]2  -  [AD  -  H  AD]T  [AD  -  H  a5] 
J=1  n_in~ 


(3-103A) 


where  from  the  standard  solution  of  least  square,  AD  is  given  by 


AD  =  (HT  H)'1  HTAD 


(3-103B) 


Using  Equations  3-103A  and  3-103B,  we  can  rewrite  Equation  3-103  as  the 


fol lowing . 


!L  RT  C 


MAX 

j=i  VcT  C 


-  K(ADT  [I  -h(ht  h)"1  HT]AD) 


(3-104) 


where  the  H  matrix  in  our  case  is: 


*  The  ir 
as  the 


initial  time  delay  difference  Adp  can  be  added  to  the  vector  A8T 
he  set  of  unknown  parameters  to  ce  estimated. 


and  I  is  the  identity  matrix. 


Also  note  that 


Ad,  =  (hT  H)'1  HT  AD 


(3-105) 


and  an  estimate  of  the  target  time  delay  difference  rate  is  given  by  : 


id,  -2d2t(  .dj 


(3-106) 


and  the  estimate  of  the  target  time  delay  acceleration  is  given  by  : 


Ad  =  2  d2  . 


(3-107) 


*  The  time  delay  difference  rate  and  acceleration  are  referenced  to 
the  measurement  averaging  time  T  seconds. 
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The  parameter  K  represents  the  functional  relationship  between  the 
ML  and  the  a  priori  information.  Unfortunately,  the  approximate  deriv¬ 
ation  of  the  ML  term  has  obscured  the  functional  form  of  K.  Simulation 
will  be  required  to  determine  the  functional  form  and  appropriate 
values  for  K. 


A  key  concept  to  understand  is  that  Equation  3-103  or  3-104  is  not 
equivalent  to  calculating  a  set  of  ML  estimates  AD  and  then  performing  a 
quadratic  LMSF  to  refine  the  estimate.  The  update-by-update  interaction 
between  the  ML  estimator  and  the  a  priori  information  yields  a  more 
accurate  estimate  and  ensures  stability  at  low  SNR. 


As  stated  earlier,  the  normalized  measured  data  vector  Rj  will 
have  a  normalization  factor  of  As/(^s  +  N$)  or  signal  power  divided  by 
signal  plus  noise  power.  Therefore,  at  low  SNR  [N$  >  >  A$] ,  the  ML 
estimate  will  be  proportional  to  SNR  (A$/Ns).  At  high  SNR  (A$  >  >  Ns), 
the  ML  estimate  will  approach  a  constant.  The  ML  estimator  can  be 
modified  to  maintain  an  approximately  linear  relationship  with  SNR 
(A$/Ns)  by  the  following  development. 


(3-108) 


•  V  V  ■ 

•S-v 


which  is  equivalent  to  the  following  by  Equation  3-37: 


The  average  SNR  estimate  for  the  block  of  data  is  given  by  Equation 
3-87  which,  after  some  algebra,  is  given  by: 

M  W/VCT  C 

M  SNRfl,;r  =  - ! - —  •  (3-110) 

AVG  *  rr~ 

m  -  w/yc'  c 

The  right-hand  side  of  Equation  3-110  is  proportional  to  the  aver¬ 
age  SNR  and  increases  as  W  increases  (0  <  W  <  M).  Therefore,  the 
right-hand  side  of  Equation  3-110  can  be  used  in  place  of  the  ML  term  in 
Equations  3-103  or  3-104  for  the  MAP  estimator. 

Wolcin6  derived  an  a  priori  knowledge  penalty  function  for  an 
assumed  piecewise  linear  frequency  trajectory  which  may  be  appropriate 
for  our  work.  His  assumptions  are  the  following: 

1.  The  densities  Pg(Adg)  and  P^(Ad^/Adg)  are  uniform  over  the 
space  of  allowable  choices  of  Adg  and  Ad^  (i.e.,  no  start  time 
delay  difference  or  time  delay  difference  rate  should  be  pre¬ 
ferred  over  any  other). 

2.  The  expected  value  of  Ad^  given  Adi  ^  and  Adi  can  be 
obtained  by  linear  extrapolation.  (The  time  delay  difference 
function  is  smooth  and  can  be  approximated  by  a  piecewise 
linear  function.) 


3.  The  sequence  {AcIq,  Ad^  ...  Aa,^}  is  second-order  Gauss-Markof f 
(note  discussion  in  Section  3.5). 


The  resulting  a  priori  knowledge  penalty  function  is  the 


fol  lowing. 


^2> 


j+2  -  2  4dj*l  *  “j1' 


(3-111] 


where  Var ( Ad^/Ad^ ,  Ad.^)  =  ac 


(3-11 1A) 


Equation  3-111a  can  be  substituted  into  Equation  3-103  or  3-T04 
for  the  a  priori  information  term. 


3.8  SUMMARY 


In  Chapter  3  we  have  derived  an  approximate  MAP  estimator 
(Equation  3-104)  for  a  target  time  difference  trajectory  as  observed 
between  two-point  sensors.  We  also  have  derived  an  SNR  estimator 
(Equation  3-37)  and  an  estimator  for  a  spectral  slope  parameter 
(Equation  3-97)  for  the  target  signal. 


Using  the  LMSF  penalty  function,  we  also  have  estimates  for  the 
target  time  delay  difference  rate  and  acceleration  (Equations  3-1C6 
ana  3-107,  respecti ve ly) . 


In  the  following  chapters  we  snail  examine  tecrniques  arc  oroolems 
in  irrplementng  the  above  estimators. 


s  V 


CHAPTER  4 


PRACTICAL  IMPLEMENTATION  OF  THE  MAP  ESTIMATORS 

4.1  INTRODUCTION 


In  Chapter  3  we  have  derived  an  approximate  MAP  estimator  for  the 
unknown  target  time  delay  difference  trajectory  due  to  a  moving  target 
over  an  extended  period  of  time  t.  Our  total  measurement  consists  of  a 
block  of  M  consecutive  measurements,  each  of  duration  T  seconds. 
(Therefore,  we  have  that  M  *  T  =  t.)  Each  of  the  M  individual  measure¬ 
ments  consists  of  the  cross-correlation  of  two-point  sensors  averaged 
for  T  seconds  and  evaluated  at  2N+1  discrete  lag  delay  times  centered 
about  zero  lag  delay. 

Equation  3-104  nominally  represents  the  MAP  estimator  for  our 
described  input  measurements.  Effectively,  the  optimum  time  delay 
difference  trajectory  is  that  which  maximizes  Equation  3-104  or  equiva¬ 
lent  over  the  block  of  input  lag  point  measurements. 

Two  general  techniques  for  solving  the  above  maximization  problem 
become  immediately  apparent.  The  first  is  to  discretize  and  limit  the 
allowable  time  delay  difference  space  IT)  (note  AD  depends  on  AST  via 
Equation  3-100).  With  only  a  finite  number  of  allowable  time  delay 
difference  trajectories,  a  solution  for  a  maximum  is  obtainable.  This 
approach  we  shall  call  the  dynamic  programming  approach. 
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The  second  approach  is  to  calculate  the  gradient  vector  of  the  MAP 


equation  (Equation  3-104  or  equivalent)  with  respect  to  the  unknown 
time  delay  difference  parameters  AST.  Using  the  gradient  vector,  a 
numerical  search  can  be  used  to  find  the  maximum.  This  approach  we 
shall  call  the  gradient  search  approach. 

In  the  subsequent  sections  we  shall  discuss  and  outline  the  above 
two  approaches  as  well  as  other  related  suboptimal  approaches. 

4.2  DYNAMIC  PROGRAMMING  APPROACH 

One  solution  to  finding  the  maximum  of  the  MAP  equation  from  a 

finite  set  of  time  delay  difference  trajectories  AD  is  to  simply  try  all  • 

of  them  and  keep  the  one  which  maximizes  the  MAP  equation.  This 

approach  is  not  practical  for  reasonable  block  sizes  and  reasonable 

7  8 

number  of  allowable  time  delay  difference  points.  Wolcin's  ’  solution 

of  a  similar  problem  for  maximizing  a  MAP  function  over  a  discretized 

space  of  complex  FFT  data  is  the  basis  for  our  approach.  Wolcin  in  his 

16 

work  used  a  dynamic  programming  technique  developed  by  Bellman. 

Wolcin's  original  algorithm  was  modified  by  LaTourette,  Greineder,  and 
g 

Wolcin  into  a  recursive  batch  technique  for  frequency  estimation  which 
parallels  our  development  for  time  delay  estimation. 
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4.2.1  Basic  Concept 


We  shall  now  briefly*  describe  the  dynamic  programming  approach 
and  pertinent  parameters  used  to  solve  our  maximization  problem  in  the 
time  delay  difference  domain. 

The  first  step  is  to  discretize  the  available  time  delay  differ¬ 
ence  space  and  limit  our  region  of  search.  The  starting  point  for  our 
discretion  process  is  the  measured  lag  point  data.  We  shall  subdivide 
the  allowable  time  delay  difference  trajectory  to  correspond  to  ' NP 1 
subdivision  between  the  discrete  lag  point  data.  The  optimum  time  delay 
difference  trajectory  shall  be  limited  to  an  initial  discretized  point 
and  one  additional  discretized  point  valid  at  the  end  of  each  of  the  M 
consecutive  measurements.  To  limit  our  search  to  a  finite  region,  we 
shall  limit  our  search  to  nominally  m  lag  points  per  update  centered 
about  our  latest  estimate.  The  m  lag  points  correspond  to  precisely 

NPT  =  m(NP  -  1)  +  1  (4-1) 


;  *  The  fine  development  of  the  dynamic  programming  algorithm  is  beyond 

f  the  scope  of  this  study.  Our  purpose  is  to  present  enough  detail  to 

!  discuss  the  recursive  modification.  For  more  information  on  this 

algorithm,  contact  the  author  or  consult  References  7,  8,  and  9. 
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discretized  points.  To  further  limit  our  search  process  we  shall 
limit  to  'NPMAX'  the  maximum  number  of  discretized  points  that  an 
optimal  time  delay  difference  trajectory  can  move  in  one  update. 

As  noted  earlier,  each  update  is  centered  about  the  latest  optimal 
estimate.  Therefore,  although  the  search  region  for  each  update 
is  precisely  NPT  points,  the  search  regions  for  adjacent  updates 
do  not  have  to  overlap.  However,  adjacent  search  regions  can  only 
be  offset  by  a  maximum  of  NPMAX  discretized  points.  This  restriction 
ensures  at  least  one  allowable  path  through  every  discretized  point 
in  the  data  block.  Figure  4-1  conceptually  illustrates  a  portion 
of  the  discretized  data  block.  Note  if  NPMAX  =  4, the  magnitude 
of  the  shift  between  update  j  +  1  and  j  +  2  is  the  maximum  allowed. 
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We  shall  now  briefly  describe  the  dynamic  programming  process. 

Let  us  assume  we  know  and  have  stored  in  matrix  t.  the  optimal  time 

J 

delay  difference  trajectory  from  each  of  NPT  discretized  points  of  time 
frame  ‘j’  back  to  our  initial  time  delay  difference  estimate  prior  to 

A 

time  frame  1.  The  matrix  T  is  defined  below. 


Tj  =  [ ADT U j ) »  A0T(2jh  ....  AdVpTj)] 
where  ADT(i.)  is  defined  as: 

J 

AD  ( i  j  )  ~  [ i g »  i^,  ^2*  •••>  i j_p  ij]  • 


(4-2) 


(4-3) 


We  shall  also  assume  that  we  have  stored  in  vector  F I T ( i  ■ )  the  sum 

J 

of  the  ML  costs  (Equation  3-84)  corresponding  to  the  optimal  trajectory 
AD(ij);  i-e., 

j 

£H(ij)  -2><V 

k=l 

where  ML(ik>  ik_^)  is  defined  from  Equation  3-83  as  the  following. 


■k-lJ 


(4-4) 


ML(tk,  1^)  = 


_CTRk 


Note  AST  =  ik  -  ik_^  . 


(4-5) 
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The  next  step  in  the  dynamic  programming  technique  is  to  find  the 
optimal  time  delay  difference  trajectory  for  all  NPT  discretized  points 
in  time  frame  j  +  1  back  to  the  original  estimate,  assuming  our  prior 
knowledge  of  the  optimal  trajectories  from  the  jth  time  frame  back  to 
the  original  estimate.  Therefore,  for  each  of  the  NPT  discretized 
points  in  time  frame  j  +  1,  we  shall  perform  the  following  calculation: 


For  each  discretized  point  1 i '  in  time  frame  j  +  1,  we  shall 


calculate  incremental  ML  functions  [ML(ij+^,  ij),  ij  :  i ^  - 
NPMAX  <_  ij  <  i^  +  NPMAX  (Equation  4-5)]  corresponding  to  the 
time  segment  from  the  end  of  time  frame  j  to  the  end  of  time 
frame  j  +  1.  From  the  above,  there  are  2  NPMAX  +  1  candidate 


optimum  time  delay  difference  trajectories  from  point  i  .+^.  These 


candidate  optimal  time  delay  difference  trajectories  shall  be 
designated  as  shown  below. 


^j+l*  ^j^  ~  M q ’  i »  ^2’  •••*  ^j»  ^j+1^  *  (^-6) 
ij  is  such  that  i^+^  -  NPMAX  <  i  .  £  i  .+^  +  NPMAX  . 

Note  that  ig,  i^,  i 2 »  •••»  ij  is  the  optimum  time  delay  difference 
trajectory  from  the  ith  point  of  the  jth  time  frame  back  to  the  initial 
time  frame. 


To  determine  the  optimal  trajectory  we  proceed  as  follows.  Each 
incremental  cost  ML(i.  ,,  i.)  is  added  to  the  accumulated  optimum  ML 
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cost  F IT ( i ■ )  corresponding  to  discretized  point  i  ■  to  form  a  candidate 

J  J 

ML  cost  F_IT(ij+i,  ij). 


'j1  ■  ML(Vr  ’j) +  • 


(4-7) 


In  addition,  using  the  candidate  optimal  time  delay  difference 


trajectories  of  Equation  4-6,  we  can  now  calculate  candidate  sum  of 


squared  residuals  SSR ( i »  ij)  about  a  quadratic  function  for  each 
candidate  trajectory.  From  Equation  3-104  we  can  write  the  following. 


SSR(ij+1,  ij)  =  ADT( i j+1 .  ij)  (I  -  H(HT  Hf1  HT)  AD(ij+1>  ij). 


(4-8) 


ij  is  such  that  i^+1  -  NPMAX  <  ij  <  i^+1  +  NPMAX 


where 


H  =  . 


(4-9) 
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The  optimal  trajectory  from  point  ij+^  can  now  be  calculated  using 
Equation  3-104  as  shown  below. 

AD( i j+1 )  =  MOj+p  ij)  for  ij  :  MAX [f_IT( i J+1 ,  ij)  -  K  SSR(ij+1,  ij)], 


iJ+1  -  NPMAX  <  ij  <  ij+1  +  NPMAX 


(4-10) 


where  K  is  a  constant  that  weighs  the  ML  function  against  the  a  priori 
information. 


This  process  is  repeated  for  all  i ^  such  that  1^  <  <  NPT^+1 

until  we  have  an  optimal  time  delay  difference  trajectory  back  from  all 
NPT  discretized  points  composing  time  frame  j  +  1. 

For  all  the  optimal  trajectories  AD  ( i  ) ,  i^  =  1,  2,  3,  ...  NPT, 
we  can  choose  the  one  with  the  largest  value  for  Equation  4-10.  This 
estimate,  i ^ ,  we  shall  use  to  center  our  next  input  lag  point  data 
measurement . 


The  dynamic  program  process  can  be  performed  inductively  until  an 
optimal  time  delay  difference  trajectory  AD ( i M)  is  determined  for  the 
entire  block  of  input  data  measurement  (M  updates).  At  that  point  we 
have  completed  the  transient  dynamic  programming  process.  From  this 
point  on  we  would  like  to  begin  the  recursive  sliding  window  update  of 
our  estimator.  By  this  we  mean  that  the  total  number  of  dat3  measure¬ 
ments  we  shall  process  after  every  update  shall  remain  fixed  at  M 


*Jk  "***  '•*  %"*  * 


updates.  As  each  new  data  measurement  is  incorporated  into  the  block  of 
data,  we  will  drop  the  oldest  data  measurement.  We  would  then  like  to 
re-estimate  the  optimal  time  delay  difference  trajectory  AD(iM+^)  for  the 
new  block  of  M  data  measurements. 


Note  that  if  we  know  the  entire  optimum  time  delay  difference 
trajectory  for  a  block  of  M  measurements,  we  effectively  have  an  esti¬ 
mate  for  the  time  delay  difference  at  any  time  inclusive  of  the  block  of 
M  data  measurements.  We  shall  nominally  choose  the  middle  time  frame 
estimate  and  the  most  recent  time  frame  estimate  for  our  outputs.  Since 
our  estimates  shall  be  recalculated  every  update,  we  shall  have  two  time 
histories  of  estimates.  The  middle  time  frame  history  shall  be  delayed 
by  half  the  block  size  from  the  most  recent  time  frame  history.  We 
would  expect  the  middle  time  frame  time  history  to  have  a  significant 
reduction  in  variance  with  respect  to  the  most  recent  time  frame 
history.  In  other  words,  we  should  expect  that  the  interpolation  pro¬ 
cess  should  be  more  accurate  than  the  extrapolation  process. 


In  the  next  section  we  shall  discuss  efficient  numerical  tech¬ 
niques  for  implementing  the  concept  of  the  recursive  batch  process. 


4.2.2  Recursive  Modification 


The  straightforward  approach  to  implement  the  recursive  batch 
process  would  be  to  simply  repeat  the  inductive  dynamic  programming 
algoritnm  over  the  entire  sliding  block  of  M  data  measurements  after 
each  update.  This  procedure  would  involve  excessive  processina  and  a 
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3.  With  Equation  4-7  in  mind,  perform  the  following  operation  to 
remove  the  effects  of  the  first  time  frame  from  FIT(i^) . 

my  =  ni(iM)  -  ML(ir  ?0)  .  (4-id 

4.  Now  simply  delete  ig  from  the  corresponding  vector  AD(iM)  in 

v 

We  are  now  in  a  position  to  add  the  effects  of  the  M  +  1  time  frame 
using  the  standard  dynamic  programming  procedure.  Candidate  optimum 
time  delay  trajectories  AD(iM+i>  y  and  their  corresponding  ML  costs 
FIT (inti*  y  can  be  obtained  via  Equations  4-6  and  4-7,  respectively. 

To  complete  our  calculation,  we  need  to  calculate  the  corresponding 
candidate  sum  of  squared  residuals  SSR(iM+p  i'M)  for  the  ft  +  1  time 
frame. 

The  straightforward  approach  to  calculate  SSR(i[vl+^,  i^)  would  be 
to  simply  re-apply  Equations  4-8  and  4-9  on  the  new  candidate  optimum 
time  delay  difference  trajectory  d>D(i..+^,  if.).  However,  Ng  ana 
Lambert^  have  derived  a  recursive  sliding  window  technique  for  updat¬ 
ing  estimates  of  polynomial  Least  Mean  Square  Fit  (LMSF )  coefficients 
which  we  can  make  use  of  to  speed  up  our  calculations.  In  their  work, 

Ng  and  Lambert^  developed  a  recursion  for  a  parameter  vector  j_n  with 
components  1  n ( i ) ,  i  -  1,  2,  3,  ...  p  +  1,  where  p  is  the  order  of  the 
polynomial  fitting  function  and  n  is  the  update  index.  In  steady  state 


the  recursion  assumes  that  a  polynomial  LMSF  is  being  performed  on  a 
sliding  data  window  of  fixed  data  length  M  +  i.  For  our  appl icatier, 


(4-12) 


have  defined  as  follows. 

Lf)+1  =  H  iM) 

where  the  matrix  H  is  defined  in  Equation  4-9  for  j  =  M. 

Let  us  assume  that  we  have  L^,  then  Ng  and  Lambert*®  give  the 
★ 

following  recursion  to  calculate  the  components  of  L^, . 


IW<1>  *  ¥‘> 

+  ^Vi  ■  V 

(4-13) 

»h<2>  ■  W2> 

(4-14) 

n.l<3>  *  V3) 

♦  <»  *  u2  V!  ♦  2  l,w(2)  -  Vld)  ' 

(4-15) 

Note  that  the  above  recursion  depends  on: 


1.  The  oldest  data  point  ig 


2.  The  newest  data  point  i^ 

3.  The  previous  values  of  L^. 


Using  Equation  4-12,  we  can  rewrite  Equation  4-3  to  the  following 


form. 


*  Our  sliding  window  has  M  data  measurements  but  M  +  1  estimates  (note 
ig  is  the  initial  estimate).  However,  we  begin  our  LMSF  at  t  =  0, 

not  t  =  1.  Equations  4-13,  4-14,  and  4-15  reflect  these  differ¬ 
ences  . 


Ln 


4-13 


SSR(i 


M+l! 


-AS.  ^M+r 


V  M<W 


•’4-16) 


T  -1 

The  matrix  (H  H)  in  our  application  is  a  <ncwn  crecalc  ul  ated 
3x3  constant  matrix.  The  vector  can  be  recurs’vely  calculated 
from  L^|  by  Equations  4-13,  4-14,  and  4-15.  The  only  quantity  remaining 
in  Equation  4-16  to  be  recursively  updated  is  i . , . , ,  i..)  — D  ( i . .  .  ,  i,,). 

This  quantity  can  be  rewritten  as  follows: 


J/T 


T/7 


—  ^M+l’  ^  — ^M+l’  W  =  ~  +  (’m+1^  "  ^ 


M+l’  M 


(4-17) 


which  completes  our  recursion  for  Equation  4-16. 

We  can  now  complete  our  process  for  the  incorporation  of  informa¬ 
tion  from  time  frame  M  +  1  by  evaluating  Equation  4-10  for  i,^  =  1, 

2,  3,  ...,  NPT. 

The  recursive  procedure  discussed  in  this  section  represents  a 
significant  reduction  in  computation  via  an  efficient  use  of  previously 
calculated  information  and  thus  makes  feasible  the  conceDt  of  a  prac¬ 
tical  recursive  batch  process. 

However,  the  strai ghtf orward  use  of  the  recursive  Datch  procedure 
may  have  potential  to  introduce  numerical  instability.  We  shall  now 
discuss  this  problem  ana  its  solution. 
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The  first  order  recursive  digital  filter  has  the  following  simple 
form  where  =1,  2,  3,  ...  N+l  is  the  input  data  sequence. 

Vl  s  “  YN  +  Y  ZN+1  •  (4-18> 

The  current  output  depends  on  all  previous  inputs  Z^ , 
i  =  0,  1,  2,  ...  N+l.  For  this  reason  these  filters  are  called 
infinite  impulse  response  filters.  If  a  <  l,*7  the  effects  of  old 
inputs  and  corresponding  errors  will  decay  exponentially,  thus  ensuring 
stability. 

Recursive  batch  filters  depend  on  a  fixed  finite  number  of  pre¬ 
vious  input  data  points.  For  this  reason  these  filters  are  called 
finite  impulse  response  filters.  These  algorithms  are  made  recursive 
at  each  update  by  removing  the  effects  of  the  oldest  data  point  and  then 
adding  the  effects  of  the  newest  data  point.  Symbolically,  these  algor¬ 
ithms  can  be  written  in  the  following  form. 

YN+1  =  YN  +  ^N+l’  ZN+1-m)  (4-19) 

where 

M  =  block  size 

ZN+i  =  newest  data 

ZN+1_M  =  oldest  data 


F(-)  =  some  recursion  function. 


Numerical  problems  can  occur  in  Equation  4-19  if  what  is  input¬ 
ted  at  the  current  update  (Z^)  and  processed  for  M  updates  is  not 
exactly  removed  at  the  M  +  1  update.  This  error  is  normally  the 
result  of  computer  roundoff  errors  in  the  recursive  computation  of 
Equation  4-19.  Note  from  Equation  4-19  that  whatever  error  we  may 
have  at  update  N  will  be  fully  carried  over  to  update  N  +  1.  It 
is  very  important  to  realize  that  this  error  has  nothing  to  do  with 
any  errors  in  the  input  data  Z^.  Z^  can  be  as  inaccurate  as  wanted. 

The  only  requirement  we  have  is  that  whatever  is  inputted  at  update 
N  and  processed  for  M  -  1  additional  updates  will  be  precisely  removed 
at  update  N  +  M. 

Equation  4-11  can  be  rewritten  into  the  form  of  Equation  4-19  if 
we  include  the  update  cycle  of  Equation  4-7;  i.e., 

P£L(  1  j+1 )  =  £11(1  j)  +  1j)  -  ML(ij+1_M,  ij^)]  .  (4-20) 

The  input  Zj+j  is  represented  by  Ml( i j+p  i j ) • 

Numerically  executing  Equation  4-20  will  result  in  computer  round¬ 
off  in  the  lower  order  bits  due  to  the  summation  process.  However, 
these  errors  should  be  small,  uniformly  distributed,  zero  mean  and 
unbiased.  Computer  simulations  have  shown  for  reasonable  computer  word 
sizes  (32  bits  on  the  VAX  11/780)  that  the  cumulative  sum  of  these 
errors  is  insignificant  for  any  reasonable  practical  application  and 
thus  can  be  ignored. 


However,  numerical  execution  of  Equations  4-13,  4-14,  and  4-15  for 

the  recursive  LMSF  parameter  results  in  the  accumulation  of  errors 

which  cannot  be  ignored.  This  fact  was  recognized  by  LaTourette  and  is 
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discussed  in  a  paper  by  Ng  and  LaTourette. 


The  recursive  equations  for  l^(j)  (Equations  4-13,  4-14, 
and  4-15)  can  be  rewritten  in  the  following  revealing  form. 


wj>  mZ  (-1)i'1  Vii 


j 

E 

i=l 

+  [(-l)j  1nJ.  j  =  1,  2,  3 


(4-21) 


where  B(^)  is  the  binomial  coefficient  defined  by 


B(  i )  =  (j  -Ji)!  i! 


(4-22) 


The  first  term  represents  integer  multiples  of  the  previous  states 
of  the  accumulation  variables  l^(j).  The  second  term  is  an  integer 
multiple  of  the  current  data  point  and  represents  the  incorporation  of 
the  current  data.  The  third  term  represents  the  removal  of  the  oldest 
data  point. 


Numerical  problems  are  encountered  when  data  incorporated  in 
term  2  are  not  precisely  removed  M  +  1  updates  later  in  term  3.  These 
errors  will  oe  accumulated  and  expanded  in  term  1.  Errors  are  primarily 
the  result  of  computer  roundoff  due  to  finite  word  sizes.  Computer 
simulations  have  indicated  the  errors  do  not  sign i f icant ly  expand  for 


the  constant  and  linear  coefficient  parameters  1M+|(1)  and  lM+^(2), 
respecti vely,  but  expand  and  become  significant  for  second  order  and 
higher  coefficient  parameters  lM+^(j),  j  >  2. 

Fortunately,  there  is  a  numerical  technique  to  overcome  this 
problem.  The  solution  is  to  simply  not  allow  any  numerical  roundoff. 
This  is  accomplished  by  representing  all  input  data  as  single  precision 
variables,  all  multiplicative  coefficients  by  integers,  and  all  accumu¬ 
lation  variables  by  double  precision  variables.  The  idea  assumes  that  a 
finite  sum  of  products  of  integers  and  single  precision  variables  will 
not  totally  use  all  the  bits  contained  in  a  double  precision  word.  Thus 
there  will  never  be  any  computer  roundoff  and  what  we  input  at  the 
current  update  will  be  precisely  removed  M  +  1  updates  later.  If  the 
local  magnitudes  of  the  input  data  i.  are  roughly  the  same  and  the  order 

J 

of  the  polynomial  LMSF  small,  the  above  assumptions  should  be  valid. 
Appendix  C  describes  an  experiment  to  demonstrate  the  above  ideas. 

Note  that  if  Equation  3-111  is  used  for  the  a  priori  penalty  func¬ 
tion  it  can  easily  be  placed  in  a  recursive  form.  Using  the  definition 
of  Equation  4-2,  Equation  3-111  can  be  written  in  the  following  form. 


=  ^-l+VV  +  ^M+1  '  2  Vi  +  Vl^  ’  (’2  '  2  V  +  V2 

(4-22) 


where 


SSDf1l i ( i f| 1 1 )  =  sum  of  the  squared  differences  of  the  optimum  tra¬ 
jectory  starting  at  point  iM+^  of  the  M  +  1  time  frame. 

Equation  4-22  involves  only  integers  and  therefore  does  not 
encounter  the  numerical  problems  previously  discussed. 

4.2.3  Initialization 

The  most  difficult  period  for  a  recursive  estimation  algorithm  is 
the  period  between  the  input  of  the  first  initial  guess  until  a  steady 
state  condition  is  reached.  This  generalization  is  certainly  true  for 
recursive  batch  algorithms  where  the  basic  assumption  is  that  you  have 
available  a  block  of  M  measurements.  The  basic  problem  is  that  many 
crucial  decisions  are  being  made  by  the  algorithm  prior  to  receiving  a 
full  complement  of  M  data  measurements. 

The  first  critical  decision  imposed  on  the  estimation  algorithm  is 
the  initial  estimate.  The  initial  estimate  may  or  may  not  be  based  on 
adequate  information.  A  poor  initial  estimate  will  start  the  algorithm 
at  an  incorrect  point  and  thus  cause  transients  and  possibly  violate 
algorithmic  assumptions.  The  above  problem  is  especially  critical,  in 
our  case,  where  we  have  assumed  a  priori  knowledge*  that  the  unknown 

*  If  Wolcin's^  a  priori  information  function  (Equation  3-111)  is  used 
this  problem  should  be  reduced.  The  sum  of  the  squared  second 
differences  is  much  more  of  a  local  process.  The  initial  point  only 
affects  the  first  squared  second  difference.  However,  an  LMSF 
quadratic  curve  fit  through  an  outlying  point  will  have  an  affect  at 
all  points. 


time  delay  difference  trajectory  is  a  quadratic  function.  An  incorrect 


initial  guess  will  appear  as  a  marked  violation  of  our  a  priori  assump¬ 


tion  and  will  remain  so  until  after  the  first  data  block  is  filled. 


Fortunately,  the  dynamic  programming  technique  used  to  implement 


our  recursive  block  estimator  has  built  within  itself  a  natural  solu¬ 


tion  to  this  problem.  This  concept  was  first  recognized  by  LaTourette 


for  this  study  and  is  described  below. 


6  7 

The  original  dynamic  programming  technique  employed  by  Wolcin  ’ 


started  at  an  initial  point  and  expanded  in  a  2NPMAX  +  1  point  tree 


structure  over  the  block  of  data  measurements.  LaTourette  introduced 


the  concept  of  limiting  the  maximum  width  of  the  tree  structure  to  NPT 


points.  Since  in  steady  state  the  basic  recursive  batch  dynamic  pro¬ 


gramming  algorithm  assumes  that  the  previous  estimate  contains  a  full 


complement  of  NPT  candidate  optimal  trajectories,  why  not  assume  NPT 


candidate  initial  estimates  centered  about  the  supplied  initial  esti¬ 


mate.  Thus,  in  essence,  we  are  estimating  the  initial  estimate  in  a 


natural  extension  of  the  basic  alaorithm. 


We  have  made  some  inroads  in  our  first  decision,  the  initial  esti¬ 


mate,  but  many  more  critical  early  decisions  must  be  made  prior  to 


filling  the  data  block  and  resulting  steady  state  operation.  For 


example,  our  a  priori  information  assumption  of  a  quadratic  optimum 


time  delay  difference  trajectory  will  have  no  affect  on  the  selection  of 


candidate  optimal  time  delay  difference  trajectories  until  the  thira 


time  frame.  Therefore,  the  first  three  points  in  all  the  cancidate 
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optimum  time  delay  difference  trajectories  and  their  corresponding 
accumulated  ML  sums  will  contain  no  assistance  from  the  a  priori 
information.  Even  the  third  time  frame  will  have  only  one  degree  of 
freedom  (4  points,  3  coefficients).  The  basic  problem  is  we  are  trying 
to  estimate  too  many  parameters  on  too  few  measurements,  which  yields 
poor  results. 


Another  view  of  this  problem  is  to  view  the  LMSF  process  as  a 

filter.  The  two-sided  bandwidth  of  the  input  ML  estimates  is  1  divided 

by  the  sampling  time  (In  our  case,  the  sampling  time  is  assumed  1  and 

therefore  the  input  two-sided  bandwidth  is  1.)  The  two-sided  bandwidth 
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of  the  LMSF  process  is  given  by  LaTourette  and  Ng  as  the  following. 


BW 


P  +  1 
N  AT 


(4-23) 


where  p  =  order  of  the  polynomial  fit 
N  =  number  of  measurements 
AT  =  sampl ing  time  =  1 
and 

BW  =  two-sided  bandwidth  of  the  LMSF  process. 


Note  if  BW  =  1,  then  the  LMSF  process  is  an  all  pass  filter. 

Examining  Equation  4-23,  it  is  apparent  that  the  higher  the  order 
‘ p '  of  the  LMSF  polynomial  fit  for  a  fixed  number  of  data  points  *  K ' , 
the  larger  bandwidth  and  therefore  less  noise  discr imination.  Con¬ 
versely,  the  lower  the  order  of  the  fit  the  more  noise  discrimination. 


However,  if  the  order  of  the  fit  is  not  adequate,  a  bias  error  will 
result  due  to  model  mismatch.  Therefore,  a  designer  must  balance  model 
mismatch  bias  errors  against  noise  discrimination.  This  trade-off 
leads  us  to  the  following  technique  for  our  initialization  problem. 

The  first  measurement  involves  the  least  amount  of  information, 
yet  we  have  two  points  to  estimate  (the  initial  point  at  the  beginning 
of  the  data  frame  and  a  second  point  at  the  end  of  the  data  frame).  For 
this  first  case  we  shall  assume  a  priori  that  our  unknown  time  delay 
trajectory  is  constant  (i.e.,  the  initial  point  is  the  same  as  the 
second  point).  For  two  points,  a  constant  (p  =  0)  should  yield  a 
minimal  model  mismatch  and  yield  our  maximum  and  only  noise  discrimin¬ 
ation  (BW  =  0.5). 

For  the  second  measurement  until  some  measurement  m  <  M^,  we  shall 
assume  that  our  optimum  time  delay  difference  trajectory  is  linear 
(p  =  1).  For  the  first  few  data  measurements  this  should  prove,  in  most 
cases,  an  excellent  model.  We  have  assumed  that  the  optimum  time  delay 
difference  trajectory  is  quadratic  over  the  entire  block  of  M  measure¬ 
ments;  therefore,  it  would  not  be  unreasonable  to  assume  that  the  opti¬ 
mum  trajectory  can  be  modeled  as  linear  for  the  first  few  time  frames. 

As  the  number  of  time  frames  ' N 1  increases,  the  noise  discriminating 
advantage  of  the  linear  fit  (p  =  1)  decreases  with  respect  to  the  qua¬ 
dratic  fit  (p  =  2),  while  the  potential  for  mismatch  bias  errors 
increases  for  the  linear  model.  The  exact  solution  for  the  optimum 
value  of  m  is  highly  dependent  on  the  actual  block  size.  A  rough  rule 
of  thumb  for  m  is  given  below. 


m  <  MIN (10,  M/2]  . 


(4-24) 


In  order  to  minimize  transients  and  smooth  the  initialization  of 
the  recursive  sum  of  squared  residuals  SSR  (Equation  4-16),  we  need  to 
update  the  matrix  (HT  H)"^-  to  correspond  to  the  number  of  data  points 
1 N '  associated  with  the  current  update.  After  M  data  updates,  the  data 
block  will  be  full  and  the  matrix  (H^  H)'^  will  therefore  become  an 
invariant  constant  matrix  for  subsequent  recursive  updates.  The 
(HT  H)"^  matrix  for  the  linear  time  delay  difference  trajectory  assump¬ 
tion  is  given  below  as  a  function  of  N.  (Note  N  equal  the  number  of  time 
frames  plus  1.) 


<hT  Linear 


u 

hl,2  " 

1,2 

h2,2  - 

where 


2(N  -  1) 
N(N  +  1) 


.  -6 

*1,2  "  N(N  +  1) 
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1 2,2  ‘  N(N  +  fj  (N  -  1) 


(4-25) 

(4-25A) 

(4-25B) 

(4-25C ) 


The  following  equations  define  the  (H1”  h)‘*  matrix  for  the  quadra¬ 


tic  time  delay  difference  trajectory  assumption. 


quadratic 


nl,2 

n2,2 

"2,3 

-hl,3 

h2, 3 

h3,3  - 

where 


9(N2  -  N)  +  6 
N(N  +  1)  (N  +  2) 


(4-26, 


_  -18(2N  -  1) 

1,2  '  N(N  +  1)  (N  +  2) 


(4-26 


,  _  30 _ 

hl,3  "  N(N  +  1)  (N  +  2) 

.  _  (2N  -  1)  (96N  -  132) _ 

n2,2  "  N(N  +  1)  (N  +  2)  (N  -  2)  (N  -  1) 


(4-26 

(4-261 


_  _ 2180 _ _ 

‘2,3  ‘  N(N  +  1)  (N  +  2)  (N  -  2) 
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_  _ 180 _ 

3,3  N(N  +1)  (N  +  2)  (N  -  2)  (N  -  1) 


(4-26 


Note  that  the  transient  recursive  updates  for  the  vector  L. 
(Equations  4-13,  4-14,  and  4-15)  are  not  affected  if  'ig1,  which 
represents  the  removal  of  old  data,  is  set  equal  to  zero.  Of  course 
only  coefficients  1(1)  and  1(2)  can  be  used  in  Equation  4-16  if  the 
linear  time  delay  difference  trajectory  is  assumed  since  (h'"  h)~^  wi 
be  a  2  x  2  matrix. 


Finally,  we  shall  complete  our  initialization  compensation  by 
applying  the  standard  correction  factor  for  the  variance  to  the 


^  ‘  -  *  \ klxl klVi  a-Vj 


a  priori  information  function  (the  sum  of  the  squared  error.)  tte  will 
accomplish  this  by  modifying  the  constant  1 K *  in  Equation  3-97  or  3-97A 
to  the  fol lowing. 


k'  =  K  (4-27) 


where 

j  =  number  of  time  frames  (note  the  number  of  data  points  is  j  +  1) 
p  =  order  of  the  polynomial  IMSF. 

Once  the  data  block  if  filled  and  our  recursive  estimator  is  in 
steady  state,  1 K '  will  be  constant  for  all  future  updates. 

4.3  GRADIENT  SEARCH  APPROACH 


An  exact  solution  for  the  parameter  vector  AD  which  maximizes  the 
MAP  equation  (Equation  3-104  or  equivalent)  over  the  block  of  M 
measurements  cannot  be  obtained  analytically.  Therefore,  the  only 
technique  to  obtain  a  precise  solution  for  aD  is  by  numerical  search. 
However,  neither  Wolcin^  nor  LaTourette^  in  their  foundation  work  in 
the  frequency  domain  had  derived  a  numerical  solution. 

4.3.1  Basic  Concept 

A  basic  tool  in  numerical  optimization  of  a  scalar  function  F  with 
respect  to  a  parameter  A  is  the  graaient  vector  defined  as: 


>1 


4-25 


fi-fr  •  (4-28) 

In  our  application  the  scalar  function  F  will  be  defined  as 
Equation  3-104*.  The  parameter  vector  A  is  by  definition  the  vector  AD 
whose  components  Ad^-  are  defined  by  Equation  3-100.  From 
Equation  3-100  we  can  define  an  alternate,  more  useful,  parameter 
vector  A  shown  below. 

AT  =  [Adg,  AS^T ,  A62T,  ...  A8mT]  .  (4-29) 

where 

Adg  =  initial  time  delay  difference  at  the  start  of  the  data  block 
A8T-  =  incremental  time  delay  difference  change  from  the  end  of 

J 

the ' j -1 '  time  frame  to  the  end  of  the  'j'  time  frame. 

3F 

Our  objective  now  shall  be  to  calculate  the  gradient  vector  from 

3A 

Equation  3-104.  To  begin  this  process,  we  shall  first  divide  the  scalar 
function  F  into  two  components. 

F (A)  =  Fj(A)  +  F2(A)  (4-30) 


*  Equation  3-104  can  be  modified  by  Equation  3-110.  However,  the 
results  obtained  in  this  section  can  be  applied  to  this  modification 
via  a  simple  application  of  the  chain  rule  for  differentiation.  The 
derivative  of  the  right-hand  expression  given  in  Equation  3-110  is 

the  following:  (M/(M-W))^  3-/3A  where  W  via  Equation  3-108  is 
the  ML  term  in  Equation  3-104,  ' M '  is  the  number  of  updates  in  the 
data  block,  and  A  is  the  unknown  parameter  vector. 


where  from  Equation  3-104  we  can  define  F^(A)  as 
M 

F:(A)  =  2-j 
j=l 

and  F2(A)  as 


(4-31 


f2(a)  =  K  ADT[I  -  H(HT  H)'1  HT]AD  . 


(4-32 


We  shall  first  calculate  the  gradient  for  Equation  4-31.  Let  us 
redefine  Equation  4-31  as  the  following. 


F i ( A )  =  f  x (A)  +  f2(A)  +...+  fM(A) 


(4-33 


where  from  Equation  3-104 


f  j  (A)  =  R.  CT(AdJ_1,  ASjT)/yCT(Adj_1,  A Q-T)  £(^dJ._1,  A&jT )  . 


(4-34 


Let  us  now  take  the  derivative  of  Equation  4-33  with  respect  to 


Adg.  Using  the  chain  rule,  we  have: 


3F,  3f  i  3f o  3 Ad  i 


M-l 

3f..  44  3  Ad, 


q  oAQq  3Ad^^i 


However,  from  our  definition  of  a  piecewise  linear  function 
(Equation  2-6),  we  have  that 


and  therefore 


3Ad . 

- L  =  1 


“"i-i 


(4-37) 


Therefore  we  can  rewrite  Equation  4-35  as  the  following. 


5F1  .  V  Sfi 

3AdQ  Z-r  3Adi_1 


(4-38) 


Now  let  us  calculate  the  derivative  of  Equation  4-33  with  respect 
to  A0iT.  Using  the  chain  rule,  we  have: 


3Fj  3fj  3f  2  9Ad1 
3A8jT  "  3A8^T  +  9Ad1  3AB]T 


(4-39) 


but  from  Equation  4-35  we  have  that 


3Ad  . 

3AB~T  =  1  ’  i  =  1,  2,  3,  ...  M 


(4-40) 


Therefore  using  Equations  4-35  and  4-40,  we  can  rewrite 


The  above  procedure  can  be  repeated  for  the  remaining  parameters 


A8i  yielding  the  following  result. 

J 

3F1  .  3fj  +  3fi 

3A8,T  3A8J  Z-r  3Ad,  , 

J  J  i=j+l  1-1 

Note  that 


(4-42) 


3F,  3f m 
I  _  m 

3A8mT  "  3A8mT  ' 


(4-43) 


Examining  Equations  4-38  and  4-42,  it  becomes  readily  apparent 
that  we  need  only  calculate  the  following  derivatives. 
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j  =  1,  2,  3,  ...  M 


(4-44) 


and 


JLl 


j  »  1,  2,  3,  ...  M 


(4-45) 


However,  these  2  •  M  derivatives  only  have  two  symbolic  forms 
nominally  specified  below  for  any  time  frame  j. 


3f 

3Adg 


and 


3f 

3A6T 


(4-46) 


Using  a  nominal  symbolic  form  for  Equation  4-34,  we  can  begin  to 
calculate  these  two  derivatives.  (The  dependence  on  the  jth  time  frame 
will  be  understood  and  not  written.) 


(C1  R) 


(4-47) 


(vn)(.'  4)  - 

yc  c 


(4-48) 


The  components  of  the  vector  £  are  defined  by  Equation  3-68  and  the 
components  of  the  vector  £  are  the  nominal  measured  lag  point  data. 
Therefore,  the  only  unknown  quantities  in  Equations  4-47  and  4-48  are 
3c1/3AdQ  and  a^/SABT  ,i  *  -N,  -N  +  1,  ...  0,  ...  N,  the  components  of 
vectors  3C/3Adg  and  3C/3A0T,  respectively.  These  quantities  we  shall 
obtain  by  differentiating  Equations  3-68  with  respect  to  Adg  and  AST, 
respectively.  Using  the  Leibniz  rule  for  differentiating  an  integral 
along  with  trigonometry  and  algebra,  we  obtain  the  following 
expressions. 

Bw  Bw 

3^.  cos(u)c  8  .j )  SincJ^-  9.)  -  cos(ajc  Adg)  Sinc(—  Adg) 


Rw 

3ci  cos(u)c  ei)  Sinc(y  -  c. 
3A8  T  =  AST 


(4-49) 


(4-50) 


where 


9i  =  T.  +  Adg  +  AST 


(4-51) 
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and  wc,  Bw,  and  T.  have  been  previously  defined  for  Equation  3-67  and 
repeated  below  for  convenience. 


0)2  +  uj 
2 


(4-52 


Bw  =  u>2  -  o)^ 


T..  =  measured  correlation  data  lag  delays. 


(4-53 


We  will  completely  specify  Equations  4-49  and  4-50  if  we  give 
their  limits  when:  AST  -*•  0;  T.  +  Adg  +  AST  ■+  0;  and  AST  -*■  0  and 


T.  +  AdQ  -  0. 


£n(^' 


«*>c  sin(o)c  ei)  Sine (-^  e.) 


«w  Qw  ** 

cos(oic  e.)  [sin(—  ei) 


Is;,  cosifs;,)! 


Bw  -2 

r  9i 


(4-54) 


aS,(*Mr)*  2 


(4-55) 


abt-o 


where 


9i  =  Ti  +  Ad0 


(4-56) 


(4-57) 


r  -  .  ■ 
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Lim 

T.+Adn+ABT-*-0 
1  U  ASTj'O 


1  -  cos(uc  s i )  Sinc(|^  ei) 
ABT 


Lim 

T  .+Adn+Af3T-*-0 
1  u  A6TV0 


(4-58) 


Lim 
A8T-0 
T  .j+Adg+0 


(4-59) 


.  $0  &  °  •  14-60) 
T^Adg-0 

We  have  now  completely  defined  the  gradient  for  Equation  4-31.  We 
shall  now  calculate  the  gradient  for  Equation  4-32.  Fortunately,  this 
gradient  can  be  easily  obtained  by  vector  calculus.  The  derivative  of 
Equation  4-32  with  respect  to  the  vector  A  can  be  written  as: 


3F2(A)  /3F2(A)^T  3AD 

3A  \  3  AO  /  3A~ 


where 

/ 3F2 (A) \  j  TIT 

\~1a5~/  =  2  K  l1  -  H(H  H)  H' ]  •  (4-62) 

Note  that  the  matrix  (I  -  H(HT  H)*^  HT)  is  symmetric  and  3AD/5A 
from  Equations  4-36,  4-37,  and  4-40  can  be  written  as  the  following 
M+l,  M+l  lower  triangular  unity  matrix 


Equation  4-61  can  now  be  rewritten  in  the  following  form. 


3^2  T 

sz~ '  &  « 

where  M  is  a  constant  M+l,  M+l  matrix  defined  below. 

T  -IT  d&P 

M  =  2  K (I  -  H(h‘  H)  H1  )  -££- 


(4-64) 


(4-65) 


and  AD  is  defined  from  the  unknown  parameter  vector  A  by  Equation  3-100. 


We  have  now  completely  defined  the  gradient  vector  of  the  MAP 

Equation  3-104  with  respect  to  the  unknown  parameter  vector  A.  With  the 

knowledge  of  the  gradient  vector,  there  are  numerous  gradient,  conju- 
21  22 

gate  direction,  and  variable  metric  algorithms  to  solve  for  the 
unconstrained  maximum. 

In  this  study,  computation  time  was  considered  critical  for  a 
practical  algorithm.  Therefore,  the  computational  method  choosen  was  a 
single  line  search  along  the  gradient  direction.  However,  in  the 


(4-68)' 


Ad^  =  Adg  +  Ag^T 

A8jT  =  ABjT,  j  =  2,  3,  4,  ...  M 
and  for  simplicity 


*WiT* 


ABmt 


(4-69) 


(4-70) 


With  the  initial  estimate  vector  A  we  now  can  obtain,  via  a  single 
line  search  along  the  gradient  direction,  the  estimate  of  the  new  opti¬ 
mum  vector  A. 


Our  algorithm  can  continue  along  the  above  lines  recursively  for 
all  subsequent  data  measurement  updates.  Effectively,  we  use  the  solu¬ 
tion  of  the  previous  update  to  yield  excellent  approximations  for  the 
current  update. 

Examining  the  above  algorithm,  it  is  apparent  that  a  given  para¬ 
meter  A8 -T  is  involved  in  M  gradient  searches  from  its  inclusion  in  the 

J 

unknown  parameter  vector  until  M  updates  later  when  it  is  used  to 

estimate  the  initial  time  delay  difference.  This  observation  leads  to 
23 

the  idea  of  modifying  our  simple  gradient  line  search  algorithm  to  a 
modified  conjugate  direction  or  variable  metric  routine  where  we  update 
the  conjugate  direction  or  metric  once  every  update  cycle.  This  repre¬ 
sents  a  significant  computational  savings  from  multiple  conjugate 
direction  updates  every  cycle.  However,  due  to  adequate  performance  of 

*  For  the  first  couple  of  updates  after  initialization,  the  previous 
time  delay  difference  change  estimate  may  not  be  our  best  choice. 

The  obvious  problem  is  an  initialization  error.  A  better  choice  for 
the  first  few  updates  may  be  AB^T  =  0. 
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the  simple  gradient  search  algorithm  and  increased  computational  costs, 
this  idea  has  not  been  pursued  for  this  study. 


Note  that  in  performing  a  line  search,  the  MAP  function 
(Equation  3-104)  must  be  evaluated  several  times  for  each  recursive 
update.  Since  the  entire  optimum  time  delay  difference  trajectory  is 
altered  each  update,  we  cannot  recursively  update  the  MAP  equation  as  we 
did  for  the  dynamic  programming  algorithm.  We  must  evaluate 
Equation  3-104  at  each  new  update. 

In  this  section  we  have  assumed  that  we  started  with  an  optimal 
estimate  after  M  time  frames.  It  is  critical  that  we  examine  how  we 
propagate  an  optimal  estimate  from  time  frame  1  through  time  frame  M. 

4.3.3  Initialization 


Note  that  after  the  first  data  measurement  the  numerical  gradient 
algorithm  has  two  unknown  parameters  to  estimate.  They  are  the  initial 
time  delay  difference  Adg  and  the  time  delay  difference  after  the  first 
time  frame  Ad^  or,  equivalently,  the  incremental  time  delay  difference 
ABjT  ( Ad =  Adg  +  A8^T).  Our  initial  estimate  for  Adg  is  the  initial 
time  delay  difference  supplied  to  the  algorithm  and  our  initial  esti¬ 
mate  for  AS^T  is  zero,  the  mean  of  all  possible  incremental  time  delays. 
Therefore,  we  shall  effectively  have  estimated  our  initial  estimate 
after  the  first  time  frame  and  shall  continue  refining  the  initial  esti¬ 
mate  after  each  subsequent  data  update  until  the  data  block  is  filled 


with  M  measurements.  The  above  refinement  process  for  the  initial  esti¬ 
mate  is  the  counterpart  of  the  multiple  hypothesis  algorithm  for  the 
initial  estimate  used  by  the  dynamic  programming  algorithm. 

Using  a  similar  argument  as  employed  for  the  dynamic  programming 
algorithm  (Section  4.2.3),  we  can  conclude  that  a  quadratic  time  delay 
difference  trajectory  model  makes  little  sense  for  the  first  few  time 
frames.  Therefore,  we  shall  assume  a  linear  time  delay  difference  model 
for  the  first  m  <  Md  time  frames  and  a  quadratic  model  afterwards.  The 
effect  of  switching  time  delay  difference  trajectory  models  and 
increasing  the  number  of  data  points  from  2  through  M  alters  the  H 
matrix.  Therefore  the  matrix  [I  -  H)”^  HT]  must  be  recalculated 

each  update  for  the  MAP  Equation  3-104  and  for  the  gradient 
Equation  4-65  until  the  data  block  is  filled  with  M  data  measurements. 
This  calculation  is  straightforward  once  the  matrix  (hT  h)'1  is  cal¬ 
culated  via  Equations  4-25  through  4-25C  or  Equations  4-26  through 
4-26F. 

Similar  to  the  dynamic  programming  algorithm,  we  can  apply  the 
standard  correction  factor  for  the  variance  to  the  a  priori  information 
function  (the  second  term  in  Equation  3-104).  We  accomplish  this  by 
using  Equation  4-27. 

Simulation  studies  involving  initializing  the  gradient  search 
algorithm,  discussed  in  Chapter  6,  indicate  additional  steps  required 
to  obtain  good  initialization  characteristics  for  the  gradient  search 
algorithm.  The  discussion  of  the  additional  initial ization  procedure 
is  contained  in  Chapter  6. 


4.3.4  Sum  of  Squared  Second  Differences  Penalty  Function 


For  completeness  and  information  we  shall  define  the  gradient 
vector  if  Equation  3-111  was  used  as  the  a  priori  informa¬ 

tion  penalty  function  instead  of  Equation  4-32. 


3F2(A) 

3AdQ  "  ° 

3F?(A) 

aSSjT-  =  ‘2K[Ad2  ‘  2  Adl  +  Ad03 
3F?(A) 

a ST"  =  -2ktAdj+i  -  3  Adj  +  3  Adj-i 
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3F  (A) 

3A8^r  =  -2K[AVl  -  2  AdM  +  *M-1> 


(4-71) 

(4-72) 

Ad.  2),  j  =  2,  3,  4,  ...  M 

(4-73) 

(4-74) 


where 

Ad-  is  defined  by  Equation  4-36 

J 

A  is  defined  by  Equation  4-29 
K  is  constant  of  proportionality. 


4.4  RELATED  RECURSIVE  BATCH  ALGORITHMS 


In  our  derivation  of  the  MAP  estimator,  several  alternate  esti¬ 


mators  became  apparent  which  are  of  theoretical  or  practical  interest. 
In  this  section  we  shall  very  briefly  identify  these  algorithms  for  the 
benefit  of  later  comparison  testing. 


The  first  alternative  algorithm  uses  Equation  3-111  for  an 
a  priori  information  function.  This  algorithm  has  been  discussed  in 
some  detail  in  Chapter  3.  We  shall  call  this  algorithm  the  Sum  of 
Squared  Differences  (SSD)  algorithms.  The  primary  purpose  of  examining 
the  SSD  algorithm  is  to  study  an  alternative  a  priori  information  func¬ 
tion.  The  SSD  algorithm  was  implemented  using  a  dynamic  programming 
approach. 

4.4.1  Sequential  Algorithm 

In  Section  3.8  we  claimed  at  low  SNR  that  maximizing  the  ML  expres 
sion  and  the  a  priori  information  expression  concurrently 
(Equation  3-104)  is  superior  to  maximizing  the  ML  expression 
(Equation  3-34)  with  respect  to  the  unknown  time  delay  difference  tra¬ 
jectory  and  then  performing  a  quadratic  LMSF  to  the  ML  estimates  to 
obtain  improved  optimal  estimates.  The  'sequential'  algorithm  performs 
the  latter  operation  using  a  gradient  search  scheme. 


4.4.2  Spectral  Estimation  Algorithms 

In  Section  3.7  an  ML  estimator  was  derived  (Equation  3-93)  to 
simultaneously  optimize  for  the  unknown  time  delay  difference  trajec¬ 
tory  and  a  spectral  parameter.  This  ML  estimator  was  combined  with  the 
SSR  a  priori  information  penalty  function  to  form  a  MAP  estimator. 
Algorithm  'SPECTRAL'  solves  for  the  above  MAP  estimator  using  a  dynamic 
programming  technique. 


In  addition,  either  the  primary  dynamic  programming  algorithm  or 
the  gradient  search  algorithm  can  estimate  a  spectral  parameter 
off-line  via  Equation  3-97. 

Many  more  recursive  batch  algorithms  or  combination  of  algorithms 
are  possible,  but  the  algorithms  discussed  in  Chapter  4  are  the  ones 
primarily  examined  in  this  study. 


CHAPTER  5 


OPTIMUM  PARAMETER  SELECTION 

The  objective  of  this  chapter  is  to  examine  the  key  parameters 
which  affect  our  two  principal  MAP  estimation  algorithms.  They  are 
the  Dynamic  Programming  (SSR)  algorithm  and  the  Gradient  Search  algor¬ 
ithm.  Both  algorithms  are  adequately  described  in  Chapter  4. 

★ 

The  key  parameters  we  shall  address  in  this  study  are  listed 
below: 

1.  K  The  proportionality  constant  between  the  ML  estimate 

and  the  a  priori  information. 

2.  M,  T  The  number  of  data  measurements  in  a  data  block  and 

the  integration  time  per  measurement. 

3.  N  The  number  of  measured  lag  points  calculated  for  each 

data  measurement. 

4.  Tau  The  time  delay  spacing  between  lag  points. 

*  The  dynamic  programming  algorithm  has  other  parameters  relating 
to  the  discretation  of  the  allowable  time  delay  difference  space. 

The  equivalent  parameters  were  investigated  by  LaTourette,  Greineder, 

and  Wolcin. ^ 


5.1  TEST  METHOD 


Our  investigation  into  optimum  parameter  selection  shall  be  by 
simulation.  We  shall  vary  a  key  parameter  while  holding  the  others  con 
stant.  Our  measures  of  optimality  will  in  general  be  RMS  time  delay 
difference  error  plotted  versus  the  key  parameter  or  parameters. 

In  order  to  help  ensure  independence  of  the  simulation  program 

from  the  derived  algorithms,  we  have  used  a  stimulation/simulation 
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program  developed  by  H.  F.  Jarvis,  Jr.  For  our  purposes  the  key 
aspects  of  the  stimulation/simulation  program  are  listed  below. 

1.  The  program  simulates  programmed  target  motion. 

2.  The  program  provides  the  measured  lag  point  data  at  the 
desired  lag  times. 

3.  The  SNR  can  be  set  to  any  level  and  held  constant  thoughout 
the  simulation. 

4.  The  simulation  program  models  the  target  spectrum  by  a  piece- 
wi se  1  inear  model . 

5.  The  program  assumes  oversampling  of  the  target  and  noise 
spectrums  and  thus  models  correlated  noise  components  for 
the  measured  lag  point  data.  The  noise  correlation  matrix 
N  takes  on  the  following  form. 

N(i,j)  =  B  fsNR  T  tC0RR(i-j)  +  KSNR  C0RR(i  +  j  +  N-l) ]  (5-1) 


CORR(t)  =  cos[ujc  tau  t] 


sin[-jr  tau  t] 


fa  tau  t 


(S-1A) 


KSNR  ■  (tTInr) 


(5-1B) 


o>c  =  the  center  frequency  (radians) 


Bw  =  the  bandwidth  (radians) 


Tau  =  the  lag  point  spacing 


i,j  =  the  lag  point  indices 


N  =  the  total  number  of  lag  points 


SNR  =  halfbeam  SNR  ratio  (linear) 


T  =  integration  time  per  measurement. 


The  right-hand  side  of  Equation  5-1A  is  a  correlation  function  for 


a  flat  broadband  power  spectrum.  The  noise  power  spectrum  has  the  same 


band  limits  as  the  flat  target  signal  power  spectrum.  Therefore, 


Figure  5-1  for  the  auto  correlation  function  of  the  target  signal  spec¬ 


ifies  Equation  5-1A  the  noise  correlation  function. 


For  all  simulations,  unless  otherwise  stated,  the  following 


shal 1  be  true . 


The  signal  and  noise  power  spectrum  shall  be  flat. 


The  SNR  shall  be  held  constant  throughout  the  simulation. 


The  noise  component  of  the  measured  lag  dat3  shall  be  correlated. 


The  MAP  estimation  algorithm  shall  not  prewhiten  the  measured 


lag  point  data. 
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5.  The  nominal  relationship  between  lag  point  spacing  and 
the  auto  correlation  function  resulting  from  a  stationary 
flat  target  power  spectrum  is  shown  in  Figure  5-1. 

6.  All  time  delay  difference  errors  shall  be  normalized  by  the 
lag  point  spacing  and  therefore  shall  be  in  terms  of  frac¬ 
tional  lag  point  spacing.  The  lag  point  spacing  is  shown  in 
Figure  5-1. 

7.  The  dynamic  programming  (SSR)  algorithm  will  operate  on 
a  time  delay  difference  space  which  has  been  quantized 

to  a  quarter  of  the  lag  point  spacing.  The  maximum  change 
for  the  dynamic  programming  algorithm  (SSR)  is  twice  the 
lag  point  spacing  per  measurement  update. 

•8.  The  nominal  data  block  size  M  is  10. 

9.  The  nominal  integration  time  T  for  each  data  update  is 
5  seconds. 

10.  The  RMS  time  delay  difference  error  is  for  the  estimate  at  the 
midpoint  of  the  data  block. 

11.  The  nominal  non-stationary  dynamic  target  time  delay  differ¬ 
ence  trajectory  is  shown  in  Figure  5-2. 

5.2  THE  PROPORTIONALITY  CONSTANT  K 

An  important  question  to  ask  is  what  should  be  the  relative 
proportionality  (K)  between  the  ML  estimate  and  the  a  priori  informa¬ 
tion  function.  The  ML  estimate  is  roughly  proportional  to  SNR  and 
the  a  priori  information  penalty  function  is  roughly  proportional 
to  inverse  SNR. 


5-7 


A  practical  view  of  the  proportionally  constant  K  is  to  view  it  as 
a  measure  of  confidence  in  the  individual  ML  estimates.  If  K  is  large, 
more  weight  is  given  to  conforming  to  the  time  history  of  ML  measure¬ 
ments.  If  K  is  small,  more  weight  is  given  to  the  individual  ML  esti¬ 
mate.  To  a  certain  extent  the  ML  estimate  and  the  a  priori  information 
scale  themselves.  At  high  SNR  the  ML  estimate  will  dominate,  which  is 
proper,  and  at  low  SNR  the  a  priori  information  will  play  an  ever 
increasing  role.  Therefore,  a  second  question  is  how  does  the  propor¬ 
tionality  constant  K  depend  on  SNR. 


We  shall  examine  the  above  questions  first  with  the  dynamic  pro¬ 
gramming  (SSR)  algorithm.  We  shall  hold  the  SNR  constant  at  five  dif¬ 
ferent  levels  (5  dB,  -10  dB,  -5  dB,  -15  dB,  -18  dB).  For  each  SNR  level  • 

★ 

we  will  vary  the  proportional ity  constant  K.  Figure  5-3  shows  the 
results  of  the  above  simulations. 


Viewing  Figure  5-3,  the  first  point  to  note  is  that  for  SNR  = 

-18  dB  the  dynamic  programming  (SSR)  algorithm  diverges  (does  not 
track)  for  K  <  0.25.  This  is  an  indicator  that  the  a  priori  informa¬ 
tion  is  required  for  the  dynamic  programming  (SSR)  algorithm  to  main¬ 


tain  track  on  the  target  time  delay  difference  trajectory  at  very 
low  SNR. 


*  The  absolute  value  for  K  is  algorithmic  dependent.  Therefore, 

all  numbers  used  for  K  should  be  taken  in  a  relative  sense,  esDecially 
between  algorithms. 


Mb. 


The  most  contradictory  thing  about  Figure  5-3  is  that  for  K  <  0.5 
the  RMS  time  delay  difference  error  is  lower  for  the  -5  dB  case  than  for 
the  +5  dB  case.  Not  shown  in  Figure  5-3  is  that  for  K  =  0  (ML  estimate 
only),  the  RMS  time  delay  difference  error  for  the  +5  dB  case  is  0.332 
fractional  lag  point  spacings  and,  as  would  be  expected,  the  RMS  time 
delay  difference  error  for  the  -5  dB  case  is  0.425  fractional  lag  point 
spacings.  The  problem  is  that  at  high  SNR  the  RMS  time  delay  difference 
error  for  the  dynamic  programming  (SSR)  algorithm  is  dominated  by  dis- 
cretation  errors.  (The  dynamic  programming  (SSR)  algorithm  only  has 
resolution  down  to  a  quarter  of  the  lag  point  spacing.)  The  a  priori 
information  tends  to  smooth  out  the  discretation  error  for  a  moving  tar¬ 
get.  However,  since  the  a  priori  information  is  roughly  inversely  pro¬ 
portional  to  SNR,  larger  values  of  K  are  required  to  obtain  the  same 
smoothing.  Therefore,  at  higher  SNR  and  a  given  value  of  K,  the  discre¬ 
tation  error  remains  constant  while  the  effective  smoothing  is  reduced. 
For  the  gradient  search  algorithm  where  there  is  no  discretation  error 
we  would  not  expect  to  see  this. 

Note  that  i or  the  +5  dB  case  as  K  becomes  large,  the  RMS  time  delay 
difference  error  approaches  the  theoretical  limit  of  0.0722  fractional 
lag  point  spacings  for  an  approximately  linear  target  motion 
(Figure  5-2).  This  theoretical  limit  is  calculated  in  Appendix  D. 

In  general,  for  all  SNR  levels  the  RMS  time  delay  difference  error 
decreases  at  first  as  K  increases  but  eventually  increases  for  ,<  large 
(> 10 ) .  A  good  overall  choice  for  K  given  the  fixed  parameters  chooser 


is  K:  0.01  <  K  <  1.  The  above  makes  sense.  At  first,  the  addition  of 
a  priori  knowledge  decreases  the  RMS  time  delay  difference  error.  But 
in  the  limit  as  K  becomes  large,  the  a  priori  knowledge  will  totally 
dominate  the  ML  estimate.  The  a  priori  information  in  and  of  itself 
contains  no  information  of  the  actual  time  delay  difference  trajectory 
without  the  ML  estimate.  Therefore,  in  the  limit  as  K  goes  to  infinity 
there  is  no  information. 

Figure  5-4  is  the  corresponding  graph  for  the  gradient  search 
algorithm.  The  SNR  is  held  constant  at  SNR  levels  of  5  dB,  -5  dB, 

-15  dB,  -18  dB,  and  -20  dB  while  the  proportionality  constant  K  is 
varied.  Figure  5-4  is  more  like  what  we  would  expect.  At  high  SNR 
(+5  dB,  -5  dB),  the  constant  K  has  little  effect  except  for  large  values 
where  we  see  signs  of  a  degradation.  Note  that  the  RMS  time  delay 
difference  errors  for  SNR  levels  of  5  dB  and  -5  dB  are  significantly 
reduced  from  the  corresponding  RMS  time  delay  difference  errors  shown 
in  Figure  5-3.  This  is  true  because  the  gradient  search  algorithm 
yields  continuous  estimates  and  therefore  does  not  encounter  discreta- 
tion  errors  at  high  SNR.  Also  note  that,  as  would  be  expected,  in 
Figure  5-4  the  RMS  time  delay  difference  errors  for  the  +5  dB  case  are 
always  lower  than  the  RMS  time  delay  difference  errors  for  the  -5  dB 
case. 

*  As  we  shall  see  later,  this  optimum  choice  is  only  valid  for  the 
choice  of  parameter  we  have  held  fixed.  In  particular,  we  shall 
see  that  as  the  data  block  size  is  increased  (M  >  10),  this 
optimum  criteria  may  no  longer  be  valid. 


From  Figure  5-4,  we  can  see  in  general  that,  as  SNR  is  lowered,  the 
a  priori  information  yields  a  greater  reduction  in  RMS  time  delay  dif¬ 
ference  error.  In  each  case  the  RMS  time  delay  difference  error  reduc¬ 
tion  appears  to  be  reduced  monotonically  by  increasing  K  until  a  criti¬ 
cal  point  where  there  is  a  rapid  degradation  in  RMS  time  delay  differ¬ 
ence  error  and  loss  of  track  with  increasing  K.  Once  again  a  suitable 
range  for  the  proportional ity  constant  K  is  K:  .1  <  K  <  2. 

5.3  THE  DATA  BLOCK  SIZE  M  AND  INTEGRATION  TIME  T 

A  key  question  to  ask  is  how  many  data  measurement  (M)  should 
we  include  in  our  data  block,  and  the  related  question,  how  long 
should  we  integrate  each  measurement. 

We  shall  examine  the  above  question  first  with  the  dynamic  pro¬ 
gramming  (SSR)  algorithm. 

The  above  questions  are  likely  to  be  complicated  by  the  target 
time  delay  difference  trajectory.  Therefore,  to  initially  minimize 
complications  we  shall  use  a  time  delay  difference  geometry  which 
is  constant  and,  by  definition  of  our  initialization  process,  right 
on  a  discretized  time  delay  difference  point  used  by  the  dynamic 
programming  (SSR)  algorithm. 

The  first  case  we  shall  look  at  is  for  the  ML  estimator  only 
(K  =  0)  with  a  integration  time  of  5  seconds  and  an  SNR  of  -15  dB. 


Figure  5-5  shows  both  the  RMS  time  delay  difference  midpoint  and  end¬ 
point  error  output  for  block  sizes  of  M  =  2,  4,  8,  10,  16,  20,  and  30. 
(Note  for  an  integration  time  of  1  second  the  algorithm  diverged  for  the 
above  scenario.) 

Examining  Figure  5-5  we  note  that  block  size  has  almost  no  effect 

on  the  RMS  time  delay  difference  error  for  block  sizes  of  8  and  longer. 

★ 

This  implies  that  the  ML  estimator  under  static  laboratory  conditions 
does  not  significantly  depend  on  moderate  block  sizes  but  obviously 
depends  on  the  integration  time  per  update.  There  is,  however,  some 
dependence  on  adjacent  states  and  thus  minimum  block  size.  This  depen¬ 
dence  is  the  result  of  the  weak  assumption  of  a  piecewise  linear  model 
and  the  limit  on  the  maximum  allowable  time  delay  difference  change  per 
measurement  update.  Evidence  of  the  dependence  is  apparent  from  the 
reduced  RMS  time  delay  difference  error  of  the  midpoint  estimate  over 
the  endpoint  estimate  (RMS  midpoint  =  0.79,  RMS  endpoint  »  1.1)  and  the 
variation  in  RMS  time  delay  difference  error  for  block  sizes  of  2  and  4. 
If  each  update  was  truly  independent,  we  would  expect  the  same  RMS  time 
delay  difference  error  for  any  point  in  the  block.  Although  in  general 
we  will  only  plot  the  midpoint  RMS  time  delay  difference  errors,  the 
endpoint  RMS  time  delay  difference  errors  for  al 1  simulations  performed 
was  always  greater  than  the  corresponding  midpoint  RMS  time  delay 
difference  errors.  Effectively,  algorithmic  decisions  are  being  made 

*  In  a  scenario  where  significant  signal  fades  are  encountered, 

this  conclusion  may  not  be  strictly  true.  The  addition  data  block 
size  will  potentially  allow  for  smoother  transition  over  signal 
fades . 


Static  Scenario,  Dynamic  Programming  (SSR) 
(T  =  5  sec,  K  =  0,  N  =  3 1  ,SNR  =  - 1 5  dB) 


at  each  update  with  a  minimum  decaying  dependence  on  surrounding 
updates.  Therefore,  the  RMS  time  delay  difference  error  of  the  ML 
estimate  at  the  midpoint  of  each  block  is  a  local  process  with  little 
correlation  with  data  measurements  more  than  one  or  two  updates 
removed.  As  we  shall  see  later,  the  slight  improvement  in  RMS  time 
delay  difference  error  for  a  block  size  of  4  is  rather  odd.  We  would 
normally  expect  an  increase  in  RMS  time  delay  difference  error. 

The  next  logical  step  is  to  repeat  the  above  experiment  using 
a  priori  information  (k  =  0.5).  For  effect  we  shall  also  lower  the  SNR 
to  -18  dB.  Figure  5-6  shows  the  resulting  RMS  time  delay  difference 
error  versus  block  size.  Since  quadractic  a  priori  information  makes 
little  sense  for  block  sizes  of  2  and  4,  we  shall  only  process  block 
sizes  of  8  or  greater. 

Figure  5-6  shows  an  entirely  different  picture  than  Figure  5-5. 

The  block  size  of  8  diverged  while  the  block  size  of  30  effectively  has 

★ 

zero  error.  In  between  we  see  a  steady  reduction  in  RMS  time  delay 
difference  error  with  increased  block  size.  As  stated  in  Chapter  3,  the 
a  priori  information  is  what  combines  the  large  observation  interval 
model  with  the  multiple  measurement  model.  In  this  case  our  assumed 
model  matches  the  actual  target  time  delay  trajectory  perfectly,  for 
any  data  length.  Therefore  we  would  expect  and  have  confirmed  in 
Figure  5-6  that  increasing  the  observation  interval  decreases  the  vari¬ 
ance  and  thus  the  RMS  error. 

*  The  zero  error  is  due  to  the  discretized  allowable  time  delay 
difference  trajectory  and  the  fact  that  the  target  signal  is 
a  constant  lined  up  with  an  allowable  discretized  point. 


me  Delay  Difference  Error  vs  Block  Size 
Scenario,  Dynamic  Programming  (SSR) 


The  above  experiments  are  very  nice  and  confirm  our  basic  con¬ 
cepts  but  are  not  very  practical  in  the  real  world.  Therefore,  we 
shall  perform  similar  experiments  using  moderate  target  motion  as 
shown  in  Figure  5-2. 

Figure  5-7  shows  the  RMS  time  delay  difference  error  versus  block 
size  for  the  ML  estimator  (K  =  0)  for  an  integration  time  of  5  seconds 
at  three  different  SNR  levels  (-5  dB,  -10  dB,  and  -13  dB).  Figure  5-8 
shows  corresponding  plots  for  an  integration  time  of  1  second  at  SNR 
levels  of  0  dB,  -5  dB,  and  -10  dB.  Figures  5-7  and  5-8  for  the  dynamic 
target  scenario  agree  with  the  corresponding  Figure  5-5  for  the  static 
scenario.  The  results  indicate  that,  for  the  ML  estimator  only  (K  =  0), 
moderate  block  sizes  have  little  effect  on  the  RMS  time  delay  difference 
error.  This  result  is  not  suprising  since  data  measurements  are  still 
loosely  coupled  by  our  assumed  piecewise  linear  model  and  therefore 
local  variance  is  a  local  issue.  For  short  block  sizes  of  4  and  2, 
there  is  a  general  increase  in  RMS  time  delay  difference  error,  except 
at  the  lowest  SNR  level  and  at  a  block  size  of  4  where  there  is  a  slight 
decrease  in  RMS  time  delay  difference  error.  This  result,  in  general, 
confirms  our  theory  on  the  loose  correlation  of  data  points  due  to  the 
piecewise  linear  assumption.  For  short  block  sizes,  increasing  the 
block  size  has  some  affect  on  the  RMS  time  delay  difference  error  while 
at  longer  block  sizes  increasing  the  block  size  has  little  affect  on  RMS 
time  delay  difference  error.  Note,  however,  that  for  the  dynamic 
scenario  (Figure  5-7)  the  dynamic  programming  (SSR)  algorithm  RMS  time 
delay  difference  error  is  significantly  higher  at  0  dB  than  the  RMS  time 
delay  difference  error  for  the  dynamic  programming  (SSR)  algorithm  at 
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-15  dB  (Figure  5-5)  using  the  static  scenario.  In  fact,  the  dynamic 
programming  (SSR)  algorithm  diverges  at  -15  dB  using  the  dynamic 
scenario.  This  demonstrates  that  results  obtained  using  the  unrealis¬ 
tic  static  scenario  will  have  little  validity  in  practice  for  the 
dynamic  programming  (SSR)  algorithm.  Finally  we  should  note  from 
Figures  5-7  and  5-8  that  the  increased  integration  time  from  1  second 
to  5  seconds  reduces  the  corresponding  RMS  time  delay  difference 
errors.  Since  there  is  relatively  small  motion  in  5  seconds  there 
should  be  little  coherency  loss.  Therefore,  as  would  be  expected, 
longer  integration,  time  will  reduce  the  estimator  variance. 

Our  next  step  is  to  repeat  the  above  dynamic  scenario  simulations 
with  the  addition  of  the  a  priori  information.  Figure  5-9  shows  the  RMS 
time  delay  difference  error  versus  block  size  with  K  =  0.5  and  an  inte¬ 
gration  time  of  5  seconds  for  five  different  SNR  levels  (+5  dB,  -5  dB, 
-10  dB,  -15  dB,  and  -18  dB).  Our  first  observation  is  that  for  high  SNR 
(+5  dB  and  -5  dB)  the  block  size  has  little  effect.  However,  there  is  a 
significant  RMS  time  delay  difference  error  reduction  from  the  corres¬ 
ponding  RMS  time  delay  difference  error  plots  (Figure  5-7)  for  the  ML 
estimate  only  (K  =  0)  simulations.  Although  the  a  priori  information 
does  reduce  the  discretation  errors  at  high  SNR,  the  block  size  has 
little  effect  on  this  reduction.  At  low  SNR  levels,  Figure  5-9  indi¬ 
cates  that  at  first  longer  block  sizes  decrease  the  RMS  time  delay  dif¬ 
ference  error.  However,  at  the  longest  block  sizes  of  20  and  30,  the 
RMS  error  is  increased.  One  possibility  is  that  the  proportional ity 
constant  K  is  a  function  of  the  block  size.  Although  our  algorithms  are 
normalized  for  the  block  size,  anything  which  would  affect  the  validity 


or  precision  of  either  the  ML  estimates  or  the  a  priori  information 
function  will  affect  the  proportionality  constant  K.  Certainly  for 
long  block  sizes  modeling  mismatches  may  compromise  the  a  priori  infor¬ 
mation  function.  To  test  this  idea  we  shall  plot,  similar  to 
Figure  5-3,  RMS  time  delay  difference  error  versus  proportionality  con¬ 
stant  K  for  SNR  levels  of  5  dB,  -5  dB,  -10  dB,  -15  dB,  and  -18  dB.  How¬ 
ever,  in  Figure  5-10,  we  shall  fix  the  block  size  at  20  and  in 
Figure  5-11  we  shall  fix  the  block  size  at  30. 

Examining  Figures  5-3,  5-10,  and  5-11,  we  can  see  that  the  param¬ 
eter  K  is  a  function  of  block  size.  In  order  to  highlight  this  depen¬ 
dency,  we  shall  plot  together  in  Figure  5-12  the  RMS  time  delay  differ¬ 
ence  error  versus  parameter  K  at  an  SNR  =  -15  dB  for  block  sizes  of  10, 
20,  and  30.  Figure  5-12  indicates  that  as  block  size  increases  the 
optimum  value  of  K  decreases  and  RMS  time  delay  difference  error 
degrades  at  lower  values  of  K.  If  increased  block  sizes  resulted  in  a 
modeling  mismatch,  it  would  be  reasonable  to  expect  that  degradation 
due  to  over-dependence  on  the  a  priori  information  would  occur  at  lower 
value  of  K.  These  conclusions  are  not  valid  at  all  SNR  levels,  indicat¬ 
ing  that  K  is  a  complex  function  of  SNR  and  block  size. 

In  order  to  validate  the  conclusions  on  data  block  size  obtained 
from  Figure  5-9,  we  shall  now  plot  in  Figure  5-13  the  RMS  time  delay 
difference  error  versus  block  size  using  the  optimal  values  of  K  found 
in  Figures  5-3,  5-10,  and  5-11.  Although  the  increase  in  RMS  time  delay 
difference  error  for  the  block  size  of  30  is  not  as  great  at  low  SNR, 
Figure  5-13  basically  agrees  with  Figure  5-9. 


ality  Constant  K,  Dynamic  Program  SSR 
(SNR  =  -15  dB,T  =  5  sec,  N  =  31) 


Dynamic  Scenario,  Dynamic  Programming 
(SSR)  (T  =  5  sec,  Optimal  K,  N  =  31) 


GO  l 


Dynamic  bcenario,  Dynamic  Programming 
(SSR)  (T  =  1  sec.  Optimal  K,  N  =  3 1 ) 


the  RMS  time  delay  difference  error  but  a  small  but  noticeable  larger 
RMS  time  delay  difference  error  occurs  at  the  block  size  of  30.  For  the 
-15  dB  case  we  note  that  the  block  size  of  10  did  not  track  for  any 
value  of  K.  In  addition,  the  block  size  of  30  now  shows  an  RMS  time 
delay  difference  error  decrease  from  the  block  size  of  20. 

The  results  of  the  simulations  of  RMS  time  delay  difference  error 
versus  block  size  for  an  integration  time  of  1  second  are  somewhat 
inconclusive.  For  the  -15  dB  case  increased  block  size  certainly  pro¬ 
duces  convergence  in  the  estimator.  However,  at  other  higher  SNR 
levels,  increased  block  size  results  initially  in  no  change  in  the  RMS 
time  delay  difference  error  and  finally,  at  the  longest  block  size,  an 
increase  in  RMS  time  delay  difference  error. 

In  conclusion,  for  the  dynamic  programming  algorithm  tnere  are 
restrictions  on  the  choice  of  block  sizes.  It  appears  that  the  optimal 
choice  of  block  size  is  dependent  on  more  than  target  dynamics  and 
potential  model  mismatches.  Since  the  block  size  affects  the  RMS  time 
delay  difference  error  via  the  LMSF  a  priori  information  penalty  func¬ 
tion,  it  is  reasonable  to  conclude  that  numerical  problems  may  be  enter¬ 
ing  computation  of  the  a  priori  information  penalty  function  for  large 
block  sizes.  Regardless,  it  would  be  unreasonable  to  assume  that 
increasing  the  number  of  unknowns  and  calculation  in  numerical  algor¬ 
ithm  will  not  eventually  lead  to  computational  problems. 


We  shall  now  examine  the  effects  of  block  size  on  the  gradient 
search  algorithm.  Since  the  gradient  search  algorithm  yields  continu¬ 
ous  estimates,  we  will  at  least  be  relieved  of  the  complication  of  the 
discretation  error. 

Once  again  to  initially  simplify  the  experiment,  we  shall  use  the 
static  geometry  at  an  SNR  equal  to  -15  dB  and  we  shall  look  at  the  ML 
only  (K  =  0)  estimate*.  Figure  5-16  gives  the  RMS  time  delay  differ¬ 
ence  error  versus  block  size  for  the  above  experiments.  Examining 
Figure  5-16,  we  notice  that  something  is  wrong.  For  a  block  size  of  30 
there  is  a  significant  increase  in  RMS  time  delay  difference  error. 

Since  the  piecewise  linear  assumption  is  common  to  both  the  gradient 
search  and  dynamic  programming  algorithms,  we  would  expect  that  long 
block  size  should  not  effect  the  RMS  time  delay  difference  error.  How¬ 
ever,  as  will  become  more  apparent  as  we  proceed,  the  gradient  search 
algorithm  has  increasing  numerical  problems  with  block  size.  The 
obvious  problem  is  that  as  the  block  size  increases,  so  does  the  dim¬ 
ension  of  the  numerical  search.  As  expected,  we  notice  a  significant 
increase  in  RMS  time  delay  difference  error  at  very  short  block  sizes  of 
2  and  4,  and  possibly  8.  This  we  may  attribute  to  the  weak  coupling  of 
the  piecewise  linear  assumption. 

Our  next  step  is  add  a  priori  knowledge.  Therefore,  similar  to  the 
dynamic  programming  experiment,  we  shall  set  K  =  0.5  and  repeat  the 
experiment  at  an  SNR  =  -18  dB.  Figure  5-17  gives  the  RMS  time  delay 
difference  error  versus  block  size  for  the  above  experiment.  We  notice 
that  RMS  time  delay  difference  error  general ly  decreases  with  block 


size  with  a  noticeable  exception  at  a  block  size  of  30.  We  should  not 
be  deceived  by  the  apparently  improved  performance  of  the  dynamic  pro¬ 
gramming  algorithm  (Figure  5-6)  at  certain  block  sizes.  As  stated 
earlier,  the  dynamic  programming  algorithm  has  the  artificial  advantage 
of  having  the  target  always  right  on  a  discretized  point.  This  experi¬ 
ment  confirms  that  something  is  happening  to  the  gradient  search  algor¬ 
ithm  at  long  block  sizes  even  for  a  very  nice  geometry. 

We  shall  now  perform  experiments  using  the  more  realistic  dynamic 
target  geometry  shown  in  Figure  5-2.  Figure  5-18  shows  the  RMS  time 
delay  difference  error  versus  block  size  for  the  gradient  search  ML 
estimator  (K  ■  0)  for  an  integration  time  of  5  seconds  at  three  differ¬ 
ent  SNR  levels  (-5  dB,  -10  dB,  -15  dB).  Examining  Figure  5-18  we  c-an 
see  the  sharp  increase  in  RMS  time  delay  difference  error  with  block 
size.  This  sharp  increase  occurs  at  lower  block  sizes  at  higher  SNR 
levels.  A  reasonable  cause  for  the  SNR  dependency  can  be  seen  from 

★ 

taking  the  derivative  of  the  right-hand  expression  of  Equation  3-10 
(ML  term)  with  respect  to  the  parameter  vector  A. 

|^[ML  term]  =  (M/(M  -  W))2  (5-2) 


*  For  both  the  gradient  search  and  dynamic  programming  algorithms,  the 
right-hand  expression  of  Equation  3-10  is  used  as  the  ML  term  in 
Equation  3-104. 


W  =  the  ML  expression  given  in  Equation  3-108 

M  =  the  block  size 

A  *  a  parameter  vector  representing  the  target  time  delay 
difference  trajectory. 

At  high  SNR  W  approaches  M  and  therefore  Equation  5-1,  represent¬ 
ing  a  portion  of  the  gradient  vector,  becomes  very  sensitive  to  numeri¬ 
cal  problems. 

Note  in  Figure  5-18,  that  for  block  size  equal  to  20,  the  RMS  time 
delay  difference  error  actually  decreases  as  SNR  is  decreased!  A 
reasonable  cause  for  this  would  be  if  the  numerical  problems  dominated 
the  nominal  RMS  time  delay  difference  error  and  if  the  numerical 
problems  decreases  with  SNR  level. 

At  all  SNR  levels  for  short  block  sizes,  we  notice  an  increasing 
RMS  time  delay  difference  error  with  decreasing  block  size.  This  we 
attribute  to  the  weak  piecewise  linear  assumption  and  the  resulting 
loose  correlation  of  data  points. 

We  shall  introduce  a  priori  information  (K  =  0.5)  into  the  dynamic 
scenario  simulation.  Figure  5-19  shows  the  RMS  time  delay  difference 
error  versus  block  size  for  the  gradient  search  algorithm  at  SNR  levels 
of  5  d8,  -5  dB,  -10  dB,  -15  dB,  -18  dB,  and  -20  dB.  The  integration 
time  is  5  seconds. 


Search,  Dynamic  Scenario  (T  =  5  sec,  K  =  0.5,  N  =  9) 


Once  again,  examining  Figure  5-19,  we  can  see,  in  general,  an  ini¬ 
tial  decrease  in  RMS  time  delay  difference  error  with  increased  block 
size,  but  at  the  longest  block  sizes  the  RMS  time  delay  difference  error 
increases.  At  the  higher  SNRs  (+5  dB,  -5  d8),  the  increase  in  RMS  time 
delay  difference  error  occurs  at  shorter  block  sizes.  These  results  are 
effectively  identical  with  those  obtained  for  K  =  0  (Figure  5-18). 
However,  as  would  be  expected,  the  addition  of  a  priori  knowledge 
reduces  the  RMS  time  delay  difference  error. 

We  shall  now  see  what  improvement  can  be  obtained  by  using  the 
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optimum  value  of  K  for  each  block  size  at  SNR  levels  of  +5  dB,  -10  dB, 
-15  dB,  and  -18  dB.  Figure  5-20  plots  the  RMS  time  delay  difference 
error  versus  block  size  for  the  optimum  choices  of  K.  The  integration 
time  is  5  seconds.  Figure  5-20  also  gives  the  optimal  values  of  K  used. 

From  Figure  5-20  we  can  see  that  using  the  optimal  value  of  K 
significantly  reduces  the  RMS  time  delay  difference  error  at  the  longer 
block  sizes.  However,  at  the  longest  block  size  of  30  the  RMS  time 
delay  difference  error  is  always  greater  than  the  shorter  block  sizes 
for  all  SNR  levels.  At  the  higher  SNR  levels  (+5  dB,  -5  dB),  there  is 
an  increase  in  RMS  time  delay  difference  error  from  a  block  size  of  10 
to  20  to  30.  At  the  lower  SNR  levels  (-15  dB,  -18  dB),  the  RMS  time 
delay  difference  error  decreases  from  a  block  size  of  10  to  a  block  size 
of  20,  but  increases  at  the  block  size  of  30.  These  results  indicate 

*  The  values  of  K  were  the  best  found. 


that  increasing  the  weighting  of  the  a  priori  information  can  signifi¬ 
cantly  reduce  the  effects  of  numerical  problems  at  longer  block  sizes 
and  high  SNR  (+5  dB,  -10  dB).  At  low  SNR  levels  (-15  dB,  -18  dB), 
decreasing  the  weight  of  the  a  priori  information  yields  better 
results.  The  increase  in  RMS  time  delay  difference  error  at  the  longest 
block  size  of  30  indicates  potential  modeling  mismatch  or  numerical 
problems. 

We  shall  now  reexamine  our  gradient  search  algorithm  performance 
as  a  function  of  block  size  with  an  integration  time  of  1  second.  At 
this  short  integration  time  our  modeling  errors  should  be  minimum.  We 
shall  use  the  realistic  scenario  described  in  Figure  5-2  for  all  simu¬ 
lations. 

Figure  5-21  shows  the  RMS  time  delay  difference  error  versus  block 
size  for  the  ML  estimator  only  (K  =  0).  The  results  of  Figure  5-21  show 
large  increases  in  RMS  time  delay  difference  with  increased  block  size, 
especially  at  high  SNR  (-5  dB).  This  result  is  the  same  as  Figure  5-18 
for  an  integration  time  of  5  seconds.  Therefore,  as  would  be  expected, 
we  still  have  a  numerical  instability  at  long  block  sizes  and  high  SNR. 
Note,  however,  the  RMS  time  delay  difference  errors  in  Figure  5-21  are 
more  than  the  corresponding  RMS  time  delay  difference  errors  in 
Figure  5-13  (at  least  for  block  sizes  not  totally  dominated  by  numer¬ 
ical  problems).  At  short  block  sizes  we  notice  the  same  increase  in  RMS 
time  delay  difference  error  with  decreasing  block  size  as  we  observed  in 
Figure  5-18  for  an  integration  time  of  5  seconds. 


Search,  Dynamic  Scena 
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We  shall  now  add  a  priori  information  to  our  simulation  (K  =  0.5). 
Figure  5-22  shows  the  RMS  time  delay  difference  error  versus  block  size 
for  our  full  MAP  estimator  (K  -  0.5). 


Examining  Figure  5-22,  we  can  still  see  a  rather  dramatic  but 
reduced  increase  in  RMS  time  delay  difference  error  with  increasing 
block  size  at  high  SNR  levels  {+5  dB,  -5  dB).  At  -10  dB  the  RMS  time 
delay  difference  error  decreases  for  increased  block  size  except  at  the 
longest  block  size  30  where  the  RMS  time  delay  difference  error 
increases.  At  -15  dB  the  gradient  search  algorithm  diverges  at  block 
sizes  of  8  and  10  and  shows  an  RMS  time  delay  difference  error  reduction 
going  from  a  block  size  of  20  to  a  block  size  of  30.  The  above  indicates 
that  the  numerical  problems  still  exist  at  high  SNR  for  the  1  second 
integration.  However,  at  lower  SNR  the  RMS  time  delay  difference  error 
tends  to  decrease  with  block  size.  Since  the  shorter  integration  time 
would  reduce  modeling  mismatches  we  might  expect  this. 

We  shall  now  replot  Figure  5-22,  but  now  using  RMS  time  delay  dif¬ 
ference  error  corresponding  to  the  optimal  value  of  K  found. 

Figure  5-23  shows  the  results  of  the  above  simulations  for  block  sizes 
of  10,  20,  and  30.  The  optimal  values  of  K  are  shown  in  brackets  below 
the  plotted  points. 

Examining  Figure  5-23,  we  see  that  increasing  block  size  now 
decreases  RMS  time  delay  difference  error  for  all  but  the  highest  SNR 
level  (+5  dB).  The  optimal  choices  of  K  generally  indicate  that  larger 
values  of  K  are  required  at  high  SNR  and  also  at  shorter  block  sizes. 


Search,  Dynamic  Scenario  (T  =  1  sec,  K  =  0.5,  N 
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The  numerical  problems  at  high  SNR  evidently  require  more  than  the 
expected  smoothing  from  the  a  priori  information. 

Figure  5-20,  using  optimal  values  of  K  for  an  integration  time  of 
5  seconds,  always  shows  an  increase  in  RMS  time  delay  difference  error 
at  the  longest  block  size.  The  corresponding  Figure  5-23  for  the  inte¬ 
gration  time  of  1  second  effectively  shows  a  decrease  in  RMS  time  delay 
difference  error  with  longer  block  size.  From  this  we  may  conclude  that 
some  of  the  increased  RMS  time  delay  difference  errors  at  a  longer  block 
size  for  the  5  second  integration  time  are  the  result  of  modeling  mis¬ 
matches.  Note  that  at  several  SNR  levels  at  the  longer  block  sizes  of 
20  and  30  the  1  second  integration  time  has  a  lower  RMS  time  delay  dif¬ 
ference  error  than  the  corresponding  5  second  integration  time. 

In  summary,  longer  block  sizes  introduce  numerical  and  modeling 
problems  for  the  gradient  search  algorithm.  At  high  SNR  the  numerical 
problems  are  more  severe,  but  can  be  reduced  by  increasing  the  weighting 
(K)  on  a  priori  information.  Therefore,  the  optimal  choice  of  block 
size,  integration  time  is  model -dependent  with  an  upper  limit  deter¬ 
mined  by  numerical  stability.  To  determine  the  optimum  selection  will 
require  experimentation. 

We  will  now  examine  the  affect  on  our  estimators  due  to  the  number 
of  lag  points  N  utilized  in  our  ML  estimator.  Since  this  parameter  does 
not  directly  affect  our  a  priori  information  function,  we  will  perform 
this  study  on  only  the  ML  portion  (K  =  0)  of  our  estimators. 
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For  this  parameter  we  can  calculate  some  special  case  analytical 
work  to  predict  our  simulation  results.  We  will  restrict  ourselves  to  a 
high  SNR  assumption  and  to  an  ML  estimate  on  a  single  measurement.  For 
this  case  we  can  only  obtain  theoretical  estimates  for  a  single  unknown 
parameter  (initial  offset  Adg  or  slope  change  AST).  Simultaneous  theo¬ 
retical  estimates  for  both  parameters  cannot  be  obtained  since  for  a 
single  measurement  their  resolution  is  ambiguous. 

We  will  therefore  restrict  our  attention  to  obtaining  theoretical 
estimates  for  two  simple  problems.  The  first  is  to  estimate  the  mean 
and  variance  of  the  estimated  initial  offset  Adg.  For  this  problem  we 
shall  assume  that  the  target  motion  is  zero  (i.e.,  AST  =  0)  and  for  con¬ 
venience  that  the  actual  target  offset  is  as  shown  below. 

Ad  0  =  0  (5-3) 

In  the  second  problem  we  will  estimate  the  mean  and  variance  of  the 
slope  parameter  estimate  AST-  We  will  again  assume  that  the  actual 
initial  otcset  Adg  =  0.  However,  we  shall  look  at  several  constant 
choices  of  the  rate  parameter  AS  and  the  observation  interval  T. 

Although  the  above  special  cases  are  very  simplified,  it  is 
reasonable  to  expect  that  results  obtained  should  apply  to  our  more  com¬ 
plicated  estimators  at  least  at  high  SNR.  By  this  we  mean  if  varying 
the  number  of  lag  points  used  significantly  affects  our  simplified 
estimators,  it  would  be  reasonable  to  assume  a  similar  effect  on  our 
related  more  complicated  estimators. 
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We  shall  examine  first  the  simpler  problem  of  calculating  the 


statistics  of  the  initial  offset  Adg  (Problem  1).  Note,  however,  that 


the  basic  mathematics  is  effectively  identical  for  both  problems. 


For  convenience  we  shall  rewrite  our  ML  estimator  below.  From 


Equation  3-83  we  have: 


=  F(AdQ) 


(5-4) 


where 


C  =  the  normalized  vector  matching  function 
R  =  the  vector  of  normalized  measured  lag  points. 


Note  that  the  parameter  of  interest  is  the  number  of  components  in 


vectors  £  and  R.  Also  note  that  we  are  maximizing  with  respect  to  Adg 


and  that  A8T  =  0. 


About  the  optimum  (Ad  =  0)  we  can  approximate  our  estimator  Ad  by  a 


standard  Newton  approximation  given  below 
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We  shall  represent  the  vector  of  normalized  measured  lag  points  R 


by  a  scaled  version  of  our  mismatch  vector  C  plus  a  noise  vector  N  (for 


convenience  the  scale  factor  will  be  1). 


R  =  C  +  N  . 


(5-8) 


The  noise  vector  N  is  zero  mean  Gaussian  with  a  covariance  matrix 


Substituting  Equation  5-8  into  Equation  5-6  and  taking  the  expec¬ 


ted  value,  the  following  result  is  easily  obtained. 


(Ad)  “o*0 


A8T=0 


(5-9) 


\mQ) 


We  will  make  the  standard  assumption  at  high  SNR  that 


AdQ=0 


(Equation  5-7)  is  insensitive  to  the  noise  vector  N  and 


A6T=0 


has  a  nominal  value  equal  to  its  noise  free  or  expected  value.  There¬ 


fore,  Equation  5-7  can  be  rewritten  as  the  following  after  substituting 


R  =  C. 
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An  expression  for  the  components  of  the  noise  correlation  matrix  N 
25 

was  obtained  from  Jarvis  and  is  given  by  Equation  5-1. 

The  components  of  the  mismatch  vector  £  can  be  obtained  from  Equa¬ 
tion  3-68.  The  expression  is  given  below. 

(Ci^Ad  =0  =  C0S(a)cTi ^  S1nc(f^  Ti)  (5-13) 

A8T=0 

Where  ti  is  the  time  delay  associated  with  each  lag  point. 

The  components  of  the  vector  3_C/3Adg  can  be  obtained  from  Equa¬ 
tion  3-68.  The  expression  is  given  below 

( 3aa^)ad0-0  =  -  ac  s(n(“c  T(>  Sinc<r  T1> 

A8T=0 
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First  since  the  data  has  been  prewhitened,  the  off-diagonal  terms 
of  the  notse  correlation  matrix  must  be  set  to  zero. 

Also  the  vectors  C  and  3_C/3Adg  undergo  the  following  transforma¬ 
tion. 

C  »  L'1  C 

and 

3C  3C 

3  Ad  g  —  3Ad  g 

-1  * 
where  L  is  determined  by  the  following  equation. 

-1  -1T  - 

N  VNORM  =  L  1  L 

and 

NORM  -  +  SN^) 

Bw  T  SNR^ 


(5-19) 


(5-20) 


(5-17) 


(5-18) 


*  In  the  actual  prewhitening  algorithm  the  matrix  L  is  calculated  for  a 
low  SNR  assumption.  This  is  done  in  order  to  have  a  prewhitening 
matrix  independent  of  SNR  which  is  valid  at  low  SNR  where  it  may 
prove  useful. 


Bw  =  bandwidth  (in  radians) 


T  =  integration  time  per  measurement 

SNR  =  halfbeam  signal-to-noise  ratio  (linear). 

The  following  figures  (5-24,  5-25,  5-26)  show  the  results  of  our 
theoretical  predictions. 

Figure  5-24  shows  the  predicted  standard  deviation  in  the  time 
delay  difference  offset  parameter  Adg  as  a  function  of  the  number  of  lag 
points  centered  about  zero  lag  delay.  The  lag  noise  components  are 
assumed  to  be  independent.  Figure  5-24  contains  results  for  both  high 
SNR  (+5  dB)  and  low  SNR  (-15  dB). 

As  would  be  expected  for  both  SNR  levels,  Figure  5-24  shows  a 
decrease  in  standard  deviation  as  the  number  of  lag  points  is  increased. 
The  reduction  is  initially  rapid  but  asymptotically  approaches  a  con¬ 
stant  value  for  large  number  of  lag  points.  Since  the  noise  is  inde¬ 
pendent  and  the  increasing  number  of  lag  points  yield  additional  but 
decaying  information  on  the  signal,  this  result  is  as  expected. 

Figure  5-25  is  identical  to  Figure  5-24  except  that  now  the  noise 
components  of  the  lag  points  are  correlated. 

Now,  however,  the  standard  deviation  of  AdQ  effectively  remains 
constant  with  minor  variations  for  all  numbers  of  lag  points.  There  is 
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Theoretical  Time  Delay  Difference  Standard  Deviation  vs 
Number  of  Lag  Points  (Correlated  Noise  with  Prewhitening) 


no  reduction  in  variance  with  increased  number  of  lag  points.  Effec¬ 
tively  the  correlated  noise  has  offset  the  advantage  of  increasing  the 
number  of  lag  points. 

Figure  5-26  contains  results  under  the  same  conditions  as  Fig¬ 
ure  5-25  but  using  the  prewhitening  assumption. 

The  prewhitening  modifications  effectively  eliminates  the  minor 
variations  in  the  standard  deviation  of  as  a  function  of  the  number 
of  lag  points.  However,  increasing  the  number  of  lag  points  does  not 
improve  our  estimate  and  the  prewhitening  modification  shows  an  insig¬ 
nificant  improvement  over  corresponding  results  in  Figure  5-25.  There¬ 
fore,  from  the  above,  prewhitening  does  not  yield  any  gain. 

We  shall  now  present  some  simulation  results  on  our  actual  estim¬ 
ators  to  compare  to  our  theoretical  results.  Although  our  theoretical 
results  assumed  zero  target  motion  and  estimated  means  and  variance,  we 
shall  use  RMS  error  as  a  measurement  and  the  modest  target  motion  shown 
in  Figure  5-2.  Both  modifications  are  more  meaningful  for  our  actual 
estimators  and  will  not  alter  our  basic  conclusions. 

All  simulations  were  conducted  for  the  ML  estimator  only  (!<  =  3). 
Each  plotted  point  represents  a  Monte  Carlo  average  of  five  simulations 
of  20  minute  duration  each.  As  per  standard  procedure,  all  simulations 
were  conducted  with  correlated  Gaussian  noise  added  to  the  measured  lag 
point  data. 
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Figure  5-27  shows  the  RMS  time  delay  difference  error  versus  the 
number  of  lag  points  using  the  gradient  search  algorithm  for  an  SNR  of 
5  dB  and  an  SNR  of  -15  dB.  The  results  agree  with  the  theoretical  and 
show  no  significant  variation  with  the  number  of  lag  points  at  either 
SNR  level. 

Figure  5-28  shows  the  RMS  time  delay  difference  error  versus  the 
number  of  lag  points  using  the  gradient  search  algorithm  with  pre¬ 
whitening  at  SNR  levels  of  5  dB  and  -15  dB.  Once  again  the  results 
agree  with  the  theory,  showing  no  significant  variation  in  RMS  time 
delay  difference  error  with  the  number  of  lag  points.  In  addition,  the 
RMS  time  delay  difference  errors  of  the  prewhitening  algorithm  are 
effectively  the  same  as  the  standard  algorithm.  This  result  indicates 
that  pre-whitening  does  not  improve  our  estimator. 

For  completeness,  Figure  5-29  shows  the  RMS  time  delay  difference 
error  versus  the  number  of  lag  points  for  the  dynamic  programming  algor¬ 
ithm.  At  low  SNR  (-15  dB),  Figure  5-29  agrees  with  Figure  5-27 
and  5-23  and  shows  no  significant  variation  in  RMS  time  delay  differ¬ 
ence  error  with  the  number  of  lag  points  used.  At  high  SNR  (5  dB)  the 
results  appear  to  be  modal.  The  RMS  time  delay  difference  errors 
measured  fall  into  two  levels.  This  variation  in  theory  is  attributed 
to  the  quantization  effect  of  the  dynamic  program  algorithm  at  high  SNR. 

We  will  now  examine  the  statistics  of  the  slope  parameter  A3T  in 
our  second  problem.  Since  we  will  be  examining  cases  where  the  actual 


slope  parameter  AST  is  not  zero.  Equation  5-5  needs  to  be  slightly  modi 
fied  to  the  following  form. 


AST  =  AST 


ast=ast 


(5-21) 


The  analytical  work  of  this  problem  follows  that  of  Problem  1  with 
the  obvious  modification.  The  results  are  summarized  below. 


The  estimator  is  unbiased  at  high  SNR. 


E [AST  -  AST]  =  0 


(5-22) 


The  variance  of  AST  is  given  below. 


Var(AST)  = 


To  complete  the  above  expression,  we  need  an  equation  for  the  com¬ 

ae 

ponents  of  the  vector  evaluated  at  AdQ  =  0  and  A3T  =  AST.  This  we 
can  obtain  from  Equation  4-50  which  is  shown  below. 


3Ci  \  cos[ojc(xi  +  AST)]  Sinc[|*(x.  +  AST)]  -  C]. 

3AST/Ad0=0  ‘  AST 


(5-24) 


ast=ast 


If  AfiT  =  0,  the  following  is  true. 


Ad0=° 

A8T=0 


.Q  is  given  by  Equation  5-14. 


A8T=0 


(5-25) 


Using  Equation  5-23  we  can  now  make  some  theoretical  predictions 
concerning  our  ML  estimator,  the  number  of  lag  points,  the  averaging 
time,  and  the  target  time  delay  rate.  Our  unit  of  measurement  will  be 
the  standard  deviation  of  our  slope  parameter  time  delay  difference 
(A0T).  (The  standard  deviation  is  the  square  root  of  the  variance.) 

a 

The  slope  parameter  AgT  is  a  measure  of  the  change  in  the  time  delay 
difference  trajectory  from  the  beginning  of  the  measurement  interval  to 
the  end  of  the  measurement  interval. 


Figure  5-30  shows  several  plots  of  standard  deviation  of  time 
delay  difference  error  versus  time  delay  difference  rate.  Each  plot 
represents  a  different  selection  for  the  number  of  lag  points.  The  SNR 
for  all  plots  in  Figure  5-30  was  set  at  5  dB ;  the  noise  components  of 
the  measured  lag  point  data  is  correlated  with  covariance  matrix  N;  the 
integration  time  T  is  set  to  5  seconds. 

Examining  Figure  5-30,  we  notice  that  for  any  selection  of  lag 
points,  the  standard  deviations  in  time  delay  difference  error  are 
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initially  equivalent  for  low  time  delay  difference  rates  and  increase 
monotonical ly  as  time  delay  difference  rate  increases.  However,  for 
smaller  number  of  lag  points,  a  rapid  increase  in  standard  deviation  in 
time  delay  difference  error  occurs  at  a  low  time  delay  difference  rate. 
Note  that  even  for  3  lag  points  a  time  delay  difference  rate  of 
0.2  fractional  lag  point  spacing  per  second  must  be  exceeded  before 
there  is  a  significant  increase  in  standard  deviation  in  time  delay  dif¬ 
ference  error. 

The  data  plotted  in  Figure  5-31  is  for  conditions  identical  to 
Figure  5-30  except  that  the  noise  components  of  the  measured  lag  point 
data  are  assumed  uncorrelated.  This  case  is  not  realistic  and  is 
presented  for  academic  interest. 

In  essence.  Figure  5-31  agrees  with  Figure  5-30.  Increasing  the 
time  delay  difference  rate  monotonical ly  increases  the  standard  devia¬ 
tion  in  time  delay  difference  error  ( AST) .  For  a  lower  number  of  lag 
points,  a  rapid  increase  in  standard  deviation  in  time  delay  difference 
error  occurs  at  a  lower  value  of  time  delay  difference  rate.  The  dif¬ 
ference  in  Figure  5-31  is  that  increasing  the  number  of  lag  points 
reduces  the  standard  deviation  in  time  delay  difference  error  even  at 
low  time  delay  difference  rates.  In  addition,  the  increase  in  standard 
deviation  of  time  delay  difference  error  is  reduced  with  respect  to 
Figure  5-30  for  all  selections  of  the  number  of  lag  points.  The  above 
effects  can  be  attributed  to  the  additional  information  obtained  due  to 
the  uncorrelated  noise  components. 


Deviation  vs  Time  Delay  Difference  Rate 
(SNR  =  5  d8.  Uncorrelated  Noise,  T  =  5  sec) 


The  data  presented  in  Figure  5-32  is  also  for  conditions  similar 
to  Figure  5-30.  However,  the  SNR  was  set  to  -15  dB  instead  of  the  5  dB 
for  Figure  5-30.  The  lag  point  noise  components  were  assumed  corre¬ 
lated  and  the  averaging  time  T  equals  5  seconds. 

The  resulting  data  shown  in  Figure  5-32  is  similar  to  the  data 
shown  in  Figure  5-30.  The  obvious  differences  are  the  larger  standard 
deviation  in  time  delay  difference  error  and  the  higher  time  delay 
difference  rates  required  before  there  is  a  significant  increase  in 
standard  deviation  in  time  delay  difference  error.  Both  observations 
can  be  largely  attributed  to  the  dominance  of  noise  at  lower  SNR. 

The  final  plotted  theoretical  data  (Figure  5-33)  was  calculated 
with  conditions  similar  to  Figure  5-30.  The  key  difference  is  that  for 
Figure  5-33  the  averaging  time  T  is  set  to  1  second  where  in  Figure  5-30 
the  averaging  time  is  5  seconds.  Once  again  the  SNR  is  set  at  5  dB  and 
the  lag  point  noise  components  are  assumed  correlated. 

Figure  5-33  demonstrates  the  obvious  affects  of  the  shorter  inte¬ 
gration  time.  For  all  selected  numbers  of  lag  points  the  standard  devi¬ 
ations  in  time  delay  difference  error  (AST),  at  the  lower  time  delay 
difference  rates,  are  higher  than  their  counterparts  in  Figure  5-30. 
However,  for  all  time  delay  difference  rates  and  number  of  lag  points 
selected,  there  is  a  significant  reduction  in  the  increase  of  standard 
deviation  in  time  delay  difference  error  with  increasing  time  delay 
difference  rate.  In  fact,  the  standard  deviation  in  time  delay  differ¬ 
ence  error  is  effectly  equivalent  for  all  number  of  1 ag  oo’nt  chosen. 


Delay  Difference  Rate  (SNR  =  5  dB,  Corre 


The  above  figures  indicate  that,  for  a  given  SNR  time  delay  differ¬ 
ence  rate  and  number  of  lag  points,  there  may  be  an  optimum  averaging 
time.  This  concept  is  developed  in  the  data  presented  in  Figures  5-34 
and  5-35.  Figure  5-34  plots  standard  deviation  of  time  delay  differ¬ 
ence  error  versus  time  delay  difference  rate  for  various  selections  of 
averaging  time  T.  The  SNR  is  set  at  5  dB;  the  number  of  lag  points  is 
set  at  9;  and  the  noise  between  lag  points  is  correlated.  The  condi¬ 
tions  for  Figure  5-35  are  the  same  as  Figure  5-34  except  that  the  SNR  is 
set  at  -15  dB. 

Two  key  points  should  be  obtained  from  Figure  5-34  and  5-35.  The 
first  is  that  for  a  given  time  delay  difference  rate  and  number  of  lag 
points  there  is  an  optimal  averaging  time  which  is  a  function  of  SNR. 

The  second  point  is  that  excessive  average  time  (T  *  10  seconds)  is 
risky  for  any  target  scenarios  featuring  significant  dynamics. 

The  following  figures  represent  simulation  data  concerning  the 
number  of  lag  points,  averaging  time,  and  the  time  delay  difference 
rate. 

Performing  realistic  statistical  simulations  on  significant 
dynamic  target  scenarios  is  almost  as  difficult  as  actual  physical 
dynamic  testing.  Significant  dynamic  target  motion  limits  meaningful 
test  time  and  the  limited  test  time  tends  to  be  dominated  by  initializa¬ 
tion  and  transient  conditions.  Therefore,  for  the  most  part,  we  make  no 
claim  that  our  simulations  confirm  our  theoretical  predictions.  (The 
maximum  target  dynamics  that  we  shall  use  is  only  approximately 
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.13  fractional  lag  points  per  second.)  However,  we  do  claim  that  our 
simulations  do  demonstrate  a  significant  role  that  averaging  time  and 
the  number  of  lag  points  play  in  a  dynamic  scenario. 

The  following  simulation  data  all  assumes  correlated  noise  between 
lag  points.  Once  again  we  shall  only  use  our  ML  estimator  (K  *  0)  and  a 
nominal  data  block  size  of  10.  In  order  to  match  our  theoretical  pre¬ 
dictions,  our  simulation  will  be  limited  to  the  continuous  gradient 
search  algorithm.  We  shall  use  the  basic  dynamic  scenario  shown  in 
Figure  5-2.  However,  we  shall  multiply  the  speed  of  the  target  by 
integer  factors,  effectively  multiplying  the  resulting  time  delay  dif¬ 
ference  rate.  Obviously,  the  higher  the  time  delay  difference  rate  the 
shorter  the  duration  of  the  simulation.  All  our  simulated  plotted  data 
represents  the  Monte  Carlo  average  of  5  simulation  runs.  For  simula¬ 
tion  run  at  4  times  the  nominal  time  delay  difference  rate,  the  simula¬ 
tions  are  5  minutes  in  duration.  For  simulation  runs  at  2  times  the 
nominal  time  delay  difference  rate,  the  simulations  are  10  minutes  in 
duration. 

Figure  5-36  shows  the  RMS  error  of  our  time  delay  difference 
estimator  versus  averaging  time  for  2  and  4  times  the  nominal  time 
delay  difference  rate.  The  number  of  lag  points  was  set  to  3  and  the 
SNR  is  5  dB.  The  data  in  Figure  5-36  is  almost  certainly  dominated  by 
transient  conditions.  However,  it  does  demonstrate  the  adverse  effects 
of  high  time  delay  difference  rates  and  long  averaging  times.  In  addi¬ 
tion,  Figure  5-36  indicates  the  potential  of  an  optimal  match  between 
averaging  time  and  target  dynamics. 
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Figure  5-37  plots  the  RMS  error  of  our  time  delay  difference 
estimator  versus  time  delay  difference  rate  for  the  number  of  lag  points 
set  at  3  and  9.  The  SNR  is  5  dB  and  the  averaging  time  is  5  seconds. 

Once  again  the  data  in  Figure  5-37  at  the  higher  time  delay  differ¬ 
ence  rates  will  be  dominated  by  transient  condition.  However, 

Figure  5-37  does  demonstrate  the  adverse  effect  of  increasing  time 
delay  difference  rate.  In  addition,  the  data  in  Figure  5-37  hints  that 
a  small  number  of  lag  points  may  adversely  affect  our  estimator  at  high 
time  delay  difference  rates. 

The  final  set  of  simulation  we  performed  at  an  SNR  of  -15  dB  and 
for  an  averaging  time  of  10  seconds.  The  purpose  of  the  simulations  was 
to  determine  the  affect  of  the  number  of  lag  points  on  the  RMS  time 
delay  difference  error  at  low  SNR.  The  time  delay  difference  rate  was 
set  at  4  times  the  nominal  rate. 

Figure  5-38  presents  the  results  of  the  above  simulations. 

Figure  5-38  indicates  a  potential  degradation  in  RMS  time  delay  differ¬ 
ence  error  if  a  small  number  of  lag  points  is  used  and  high  target  time 
delay  difference  dynamics  are  encountered.  This  result  is  in  agreement 
with  Figure  5-37. 

Despite  quantization  error  which  complicates  the  dynamic  programm¬ 
ing  algorithm,  limited  simulation  results  indicate  similar  behavior  as 
recorded  in  Figure  5-36  through  5-38  for  the  gradient  search  algorithm. 
The  dynamic  programming  algorithm  is  hard-limited  to  the  maximum  time 
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delay  difference  rate  that  it  can  follow  by  the  parameter  NPMAX.  This 
parameter  limits  the  maximum  delay  difference  change  allowed  during  a 
single  measurement  update.  If  NPMAX  is  fixed,  the  longer  the  averaging 
time  the  smaller  the  maximum  allowable  time  delay  difference  rate. 

In  summation,  we  have  determined  the  following  conclusions  con¬ 
cerning  the  number  of  lag  points  N  and  the  averaging  time  T.  Our 
conclusions  will  be  limited  to  the  realistic  assumption  of  correlated 
noise  between  measured  lag  points. 

1.  For  low  target  time  delay  difference  dynamics  and  reasonable 
averaging  times,  the  variance  of  our  ML  estimator  is 
effectively  independent  of  the  number  of  lag  points.  Since 
processing  additional  lag  points  is  computationally  costly  we 
would  want  to  minimize  the  numbers  of  lag  points  (3)  pro¬ 
cessed. 

2.  For  high  target  time  delay  difference  dynamics  and  reasonable 
averaging  times,  there  is  potential  degradation  in  our  ML  for 
a  very  small  number  of  lag  points  (3). 

3.  Combining  conclusions  1  and  2,  a  recommended  number  of  lag 
points  is  either  5,  7,  or  9. 

4.  Depending  on  the  severity  of  the  target  tine  delay  difference 
dynamics,  there  is  an  optimum  averaging  time.  Reasonable 
averaging  times  appear  to  be  in  the  range  of  1  to  5  seconds. 
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The  final  parameter  of  interest  is  the  lag  point  spacing  Tau.  We 
can  use  Equation  5-23  to  make  some  theoretical  predictions  about  the 
effects  of  lag  point  spacing  on  our  theoretical  ML  estimator.  In  the 
theoretical  predictions  and  simulations  to  follow,  we  shall  use  multi¬ 
ples  of  our  basic  lag  point  spacing  illustrated  in  Figure  5-1. 


In  our  first  theoretical  analysis,  we  shall  calculate  standard 
deviation  of  a$T  (time  delay  difference  change  over  the  measurement 
interval)  versus  lag  point  spacing  (based  on  multiples  of  the  basic  time 
delay  between  lag  points).  The  above  calculation  was  performed  for 
several  time  delay  difference  rates  varying  from  0  to  1  fractional  lag 
point  spacing  per  second.  In  all  cases,  the  fractional  lag  point 
spacing  is  referenced  to  our  nominal  lag  point  spacing.  The  SP)R  is  set 
to  5  dB,  averaging  time  T  is  5  seconds  and  the  number  of  lag  points  N  is 
set  at  9.  The  noise  between  lag  points  is  correlated. 


Figure  5-39  illustrates  the  above  calculations.  Examining 

★ 

Figure  5-39  we  notice  that  for  all  time  delay  difference  rates  the 
standard  deviation  of  our  time  delay  difference  estimate  degrades  for 
lag  point  spacing  multiples  greater  than  twice  our  nominal  lag  point 
spacing.  This  result  is  not  surprising  after  looking  at  Figure  5-1 
which  shows  the  relationship  between  the  nominal  lag  point  spacing  and 
the  nominal  target  signal  auto  correlation  function.  Sampling  at  4  or 
more  times  the  nominal  lag  point  spacing  will  yield  at  most  1  lag  point 
between  the  first  two  zero  crossings.  Therefore,  we  are  obviously 
undersampling  which  results  in  loss  of  information. 

*  Unplotted  points  represent  invalid  results  for  our  theoretical 
estimate  of  standard  deviation. 
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At  moderate  time  delay  difference  rates  we  notice  from  Figure  5-39 
that  decreasing  the  lag  point  spacing  from  twice  the  nominal  value  has 
very  little  affect  on  our  estimator.  As  time  delay  difference  rates  are 
increased,  lowering  the  lag  point  spacing  results  in  an  increasing 
degradation  of  our  time  delay  difference  estimator  standard  deviation. 

Overall  from  the  theoretical  results  presented  in  Figure  5-39  we 
can  conclude  that  either  1  or  2  times  our  nominal  lag  point  spacing  will 
be  an  optimum  choice. 

Figure  5-40  is  a  repeat  of  Figure  5-39  except  that  the  number  of 
lag  points  is  set  to  3.'  The  general  results  for  Figure  5-40  are  the 
same  as  Figure  5-39.  However,  the  noticeable  improvement  in 
Figure  5-39  at  high  time  delay  difference  rates  for  1  or  2  times  the 
nominal  lag  point  spacing  is  noticeably  reduced  in  Figure  5-40.  This 
result  agrees  with  our  earlier  results  indicating  that  at  high  time 
delay  difference  rates  too  few  lag  points  will  cause  a  degradation  in 
our  ML  estimator. 

The  final  figure  is  for  academic  purposes.  The  theoretical  pre¬ 
dictions  are  for  conditions  identical  to  Figure  5-39  except  that  the 
noise  between  lag  points  is  uncorrelated. 

In  Figure  5-41  the  standard  deviation  of  our  time  delay  difference 
estimator  is  degraded  at  small  lag  point  spacings  for  all  time  delay 
difference  rates.  There  is  a  clear  optimal  lag  point  spacing  ranging 
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Figure  5-41.  Theoretical  Time  Delay  Difference  Standard  Deviation  vs  Lag 
Point  Spacing  (SNR  =  5  dB,  Uncorrelated  Noise,  N  =  9,  T  =  5) 


from  0.5  times  the  nominal  lag  point  spacing  at  low  time  delay  differ¬ 
ence  rates  to  2  times  the  nominal  lag  point  spacing  at  high  time  delay 
difference  rates. 

Theoretical  predictions  calculated  at  low  SNRs  essentially  yield 
identical  results  as  Figures  5-39  through  5-41  except  with  proportion¬ 
ately  higher  time  delay  difference  standard  deviations.  These  results 
will  not  be  presented. 

We  will  now  attempt  to  confirm  our  theoretical  prediction  with 
simulations  on  our  major  algorithms.  Since  high  time  delay  difference 
rates  result  in  questionable  simulation  results,  we  will  only  perform 
simulations  with  our  nominal  low  dynamic  target  geometry  (Figure  5-2). 

All  simulation  results  represent  the  average  of  5  Monte  Carlo 
simulation,  each  of  20  minutes  duration.  The  averaging  time  is  set  to 
5  seconds,  the  number  of  lag  points  is  set  to  9,  and  the  noise  between 
measured  lag  points  is  correlated.  Once  again  we  shall  limit  our  simu¬ 
lations  to  our  ML  estimator  only  (K  =  0).  As  always,  all  RMS  time  delay 
difference  errors  are  referenced  to  our  nominal  lag  point  spacing. 

Figure  5-42  presents  the  results  of  the  above  simulations  for  the 
gradient  search  algorithm  at  SNR  levels  of  5  dB  and  -15  dB. 

At  +5  dB  the  RMS  error  of  our  time  delay  difference  estimator 
varies  wildly  with  lag  point  spacing.  There  is  a  distinct  optimum  at 
the  nominal  lag  point  spacing.  Fortunately,  the  -15  dB  simulation 
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results  essentially  agree  with  our  theoretical  predictions.  At  short 
lag  point  spacings  there  is  small  variation  in  the  time  delay  difference 
RMS  error.  At  lag  point  spacings  greater  than  our  nominal  lag  point 
spacing,  there  is  a  gradual  increase  in  time  delay  difference  RMS  error 
until  at  lag  point  spacing  of  8  times  our  nominal  lag  point  spacing  the 
gradient  search  algorithm  diverges.  The  erratic  results  obtained  at 
5  dB  we  will  attribute  to  numerical  problems  that  have  been  examined 
ear  1 ier. 


The  dynamic  programming  algorithm  complicates  our  simulation  due 
to  its  quantization.  Reducing  the  lag  point  spacing  automatically 
reduces  lag  point  resolution  and  the  maximum  time  delay  difference 
change  allowed  in  a  given  measurement. 


Figure  5-43  presents  simulation  results  for  the  dynamic  program¬ 
ming  algorithm  at  SNR  levels  of  5  dB  and  -13  dB.  The  conditions  are 
identical  to  the  gradient  search  algorithm  shown  in  Figure  5-42. 


The  results  shown  in  Figure  5-43  are  interesting.  At  the  lowest 
lag  point  spacing,  the  maximum  time  delay  difference  update  (NPMAX)  is 
about  the  order  of  the  actual  time  delay  difference  rate.  This  should 
result  in  a  degradation  in  RMS  time  delay  difference  error  as  seen  in 
Figure  5-43.  At  the  high  undersampled  lag  point  spacing  the  RMS  time 
delay  difference  error  increases  and  the  algorithm  diverges  as  would  be 
expected.  For  moderate  lag  point  spacings,  Figure  5-43  shows  the  RMS 
time  delay  difference  errors  for  the  5  dB  case  remaining  constant 
whereas  an  increase  in  RMS  time  delay  difference  error  c3n  be  seen  from 
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1/4  the  nominal  lag  point  spacings  to  2  times  the  nominal  lag  point 
spacing.  The  increase  in  RMS  time  delay  difference  error  with. increas¬ 
ing  lag  point  spacing  would  be  reasonable  since  the  resolution  of  the 
time  delay  differences  is  decreasing  linearly  with  increasing  lag  point 
spacing.  At  high  SNR,  numerical  resolution  problems  may  account  for  the 
relative  constant  RMS  time  delay  difference  error  with  lag  point 
spacing. 
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One  should  not  be  deceived  by  apparent  improvement  in  the  dynamic 
programming  algorithm  performance  at  1/4  the  lag  point  spacing.  The 
dynamic  following  ability  for  the  same  complexity  algorithm  has  been 
severely  limited  by  the  resulting  smaller  maximum  allowable  time  delay 
difference  change  (NPMAX).  Increasing  NPMAX  will  result  in  signifi¬ 
cantly  more  computational  requirements.  In  the  final  analysis,  the 
a  priori  information  will  almost  certainly  reduce  any  apparent  gain  to 
zero.  Therefore,  even  for  the  dynamic  programming  algorithm  we  will 
recommend  that  the  nominal  lag  point  spacing  be  adopted. 
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CHAPTER  6 


COMPARATIVE  ALGORITHM  PERFORMANCE 


6.1  INTRODUCTION 

In  this  chapter  we  have  two  basic  objectives.  The  first  is  to 
present  the  performance  of  our  best  MAP  estimation  algorithm  referenced 
to  a  basic  standard.  The  best  MAP  estimation  algorithm  is  the  gradient 
search  algorithm  and  the  basic  standard  was  choosen  as  the  Cramer  Rao 
Lower  Bound  (CRLB)  for  coherent  integration  on  a  stationary  signal 
(no  a  priori  information).  The  second  objective  is  to  present  relative 
performance  for  other  alternatives/approaches  discussed  in  Chapters  3 
and  4. 

The  above  objectives  were  accomplished  by  conducting  Monte  Carlo 
simulations  on  the  gradient  search  and  other  related  algorithms.  The 
results  of  these  simulations  are  plotted  for  visual  comparison.  Unless 
otherwise  stated,  all  plotted  data  points  consist  of  the  average  of 
five  20  minute  simulation  runs. 

The  simulation  data  is  presented  using  five  different  criteria  as 
outlined  below. 

1.  Semi  static  performance  of  the  gradient  search  algorithm. 

*  The  geometry  used  is  shown  in  Figure  5-2.  This  geometry  is  not 

static,  but  has  low  time  delay  difference  rates.  The  reasons  for  not 
using  a  strictly  static  geometry  are  discussed  in  Chapter  5. 


2.  Semi-static  performance  comparison  between  the  gradient  search 
algorithm  and  other  related  MAP  algorithms. 

3.  Performance  of  the  gradient  search  algorithm  with  target 
dynamics. 

4.  Performance  of  the  gradient  search  algorithm  with  spectral 
mismatch. 

5.  Initialization  performance  of  the  gradient  search  and  dynamic 
programming  algorithms. 

6.2  SEMI  STATIC  PERFORMANCE  OF  THE  GRADIENT  SEARCH  ALGORITHM 

The  time  delay  difference  trajectory  shown  in  Figure  5-2  was  used 
for  all  simulations  in  this  section. 

In  order  to  minimize  the  effects  of  initialization  on  steady  state 
algorithm  performance,  all  simulations  were  artificially  updated  with 
the  correct  trajectory  for  the  first  10  updates  (block  size).  From  the 
11th  update  onward  the  algorithm  proceeded  normally.  This  procedure 
emphasizes  a  key  point.  Tracking  and/or  estimation  algorithms  are 
capable  of  maintaining  track  of  valid  estimates  in  steady  state  at  a 
lower  SNR  than  they  are  capable  of  initializing  on.  This  points  out  the 
important  aspect  of  track  initiation  in  cases  of  target  dynamics. 
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For  our  reference  function  we  shall  use  the  CRLB  modified  to 
reflect  standard  deviation  in  time  delay  difference  error  in  units  of 
fractional  lag  point  spacing. 


The  CRLB  bound  for  the  case  of  a  single  stationary  target  as  seen 

5 

from  two  sensors  without  a  priori  information  can  be  obtained  from  Ng 
26 

or  Carter.  Their  result  is  given  below. 


VAR(t)  >  2tt  1 2T  /  yt 


(ii)  2  . 

sfcl  “  d“ 


(6-1) 


where 


/  hh 

y  + 


Gn(w) 


(6-2) 


Gs(oj)  =  the  signal  power  spectrum  as  seen  in  identical  receive 
channels 


(co)  =  the  noise  power  spectrum  as  seen  in  identical  receive 
channels 


T  =  the  observation  interval 


the  time  delay  difference  estimate. 


We  shall  assume  that  the  signal  and  noise  are  constant  spectral 
densities  as  specified  in  Equations  2-3  and  2-4.  Therefore,  Equa¬ 
tion  6-2  can  be  rewritten  to  the  following  convenient  form. 


Using  Expression  6-1  through  6-5  we  obtain 
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1  +  2SNR 
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(6-6) 


We  will  now  modify  Equation  6-6  to  reflect  time  delay  difference 
standard  deviation  in  units  of  fractional  lag  point  spacing. 


STD(t)  =  yVAR(t)/Tau 


(6-7) 


where 


STD  =  our  modified  CRLB 

Tau  =  the  nominal  lag  point  spacing. 

We  shall  use  Equation  6-7  to  calculate  the  modified  CRLB  for 
observation  time  T  =  5  and  T  =  50.  T  =  5  is  the  observation  time  for 
our  nominal  measurement  update  and  T  =  50  is  the  observation  time  for 
our  nominal  block  size  of  10. 

Figure  6-1  consists  of  three  pairs  of  plots.  The  first  pair  is  the 
midpoint  and  endpoint  RMS  time  delay  difference  errors  versus  SNR  (dB) 
for  a  5  second  (T  =  5)  measurement  update  integration  time.  The  second 
pair  is  the  midpoint  and  endpoint  RMS  time  delay  difference  errors 
versus  SNR  (dB)  for  a  2  second  (T  =  2)  measurement  update  integration 
time.  The  third  set  of  plots  is  our  modified  CRL3  for  an  integration 


time  of  T  =  5  and  T  =  50  seconds. 


The  gradient  search  algorithm  was  used  for  all  simulations  shown 
in  Figure  6-1.  The  block  size  was  fixed  at  10  (M  =  10).  The  number  of 
lag  points  was  9  (N  =  9)  and  the  a  priori  information  proportional ity 
constant  was  set  to  1.5  (K  =  1.5). 

Examining  Figure  6-1  we  can  make  the  following  observations. 

1.  The  RMS  time  delay  difference  error  for  the  endpoint  measure¬ 
ment  is  always  qreater  than  the  corresponding  midpoint 
measurement.  At  high  SNR  this  difference  is  greater  and  at  low 
SNR  the  difference  is  smaller.  In  between,  the  RMS  time  delay 
difference  error  difference  is  relatively  constant  at  about 

1.5  dB.  At  high  SNR  our  midpoint  and  endpoint  estimators 
should  be  dominated  by  the  weak  piecewise  linear  assumption. 

As  the  SNR  is  lowered,  our  midpoint  and  endpoint  estimators 
will  become  more  dominated  by  our  quadratic  time  delay  differ¬ 
ence  trajectory  assumption. 

2.  The  RMS  error  for  the  2  second  measurement  update  integration 
time  is  always  greater  than  the  corresponding  RMS  error  for 
the  5  second  measurement  update  integration  time.  Since  the 
time  delay  difference  rate  is  small,  this  result  is  as 
expected. 

3.  The  midpoint  RMS  time  delay  difference  error  for  the  5  second 


integration  time  is  always  less  than  the  corresponding 
5  second  CRL3  and  at  low  SNR  is  less  than  the  50  second  CRL3. 


This  result  is  not  a  contradiction  of  the  CRLB.  The  proper 
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calculation  of  the  CRLB  requires  an  additional  term  due  to 
a  priori  information.  We  have  neglected  this  term.  Even  at 
high  SNR  the  weak  piecewise  linear  assumption  is  a  priori 
information,  thus  we  would  expect  our  estimator  to  outperform 
the  T  =  5  CRLB  as  a  result  of  this  information.  As  the  SNR  is 
lowered,  the  stronger  quadratic  time  delay  difference  assump¬ 
tion  will  add  more  a  priori  information  to  our  MAP  estimator 
via  the  LMSF  penalty  function  (see  Equation  3-104).  This  is 
true  since  the  relative  magnitude  between  our  ML  term  and  our 
a  priori  information  term  will  favor  the  a  priori  information 
term  as  SNR  is  lowered.  This  a  priori  information  is  not 
accounted  for  in  our  CRLB  calculation. 

4.  Our  final  observation  is  that  our  RMS  time  delay  difference 
error  curves  for  our  gradient  search  algorithm  apparently 
deviate  from  the  expected  lower  trend  at  high  SNR.  This  may  be 
a  result  of  numerical  problems  at  high  SNR  which  is  discussed 
earl ier. 

We  will  now  examine  the  overall  effect  of  the  proportional ity 
constant  K  on  our  gradient  search  algorithm. 

Figure  6-2  presents  the  RMS  time  delay  difference  errors  versus 
SNR  for  three  different  choices  of  the  proportional ity  constant  K.  Once 
again  included  in  Figure  6-2  are  the  T  =  5  and  T  =  50  second  modiried 
CRLB  curves.  The  choices  for  the  proportional ity  constant  K  are 
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K  =  1.5,  K  =  0.75  and  K  =  0.  Note  that  K  =  0  is  our  ML  estimator  only. 

The  block  size  M  is  set  to  10  and  the  number  of  lag  points  N  is  9. 

Examining  Figure  6-2  we  can  make  the  following  observations. 

1.  At  low  SNR  (<  -15  dB),  the  K  =  1.5  option  has  the  smallest  RMS 
time  delay  difference  error.  This  fact  is  in  agreement  with 
data  analyzed  in  Chapter  5.  Since  good  low  SNR  performance  is 
desirable,  we  have  choosen  K  =  1.5  as  our  nominal  design. 

2.  Below  -18  dB  the  ML  estimator  (K  *  0)  diverges.  This  results 
states  that  strong  a  priori  information  is  required  for  main¬ 
taining  track  or  stability  for  the  given  test  conditions. 

3.  Neglecting  algorithm  stability  at  low  SNR,  Figure  6-2  indi¬ 
cates  a  relatively  small  variation  in  RMS  time  delay  differ¬ 
ence  error  between  K  =  0  and  K  =  1.5.  This  result  states  that 
the  a  priori  informaton  function  is  primarily  required  for 
stability  at  low  SNR. 

4.  The  modified  CRLB  would  be  roughly  equivalent  to  the  ML  only 
(K  =  0)  RMS  time  delay  difference  error  curve  if  we  set  the 
CRLB  averaging  time  to  25  seconds.  Therefore,  if  we  equate 
averaging  time  to  a  priori  information,  we  may  conclude  that 
the  piecewise  linear  assumption  is  equivalent  to  averaging  for 
five  5  second  measurement  updates.  This  result  would  be  in 
agreement  with  test  results  in  Chapter  5.  In  Chapter  5  the  ML 


only  estimator  (K  =  0)  indicated  variation  in  RMS  time  delay 
difference  error  for  block  size  of  2  and  4;  however,  there  was 
no  significant  variation  in  RMS  time  delay  difference  error 
for  block  sizes  greater  than  8. 

The  final  plot  in  this  section  will  present  the  performance  of  the 
SNR  estimate  for  the  gradient  search  algorithm.  Since  there  is  a  known 
bias  due  to  the  data  generation  program  at  high  SNR,  we  shall  present 
only  the  standard  deviation  in  the  SNR  estimate.  The  standard  deviation 
was  calculated  from  the  SNR  estimate  in  units  of  dB. 

Figure  6-3  presents  the  standard  deviation  of  the  SNR  estimate 
versus  SNR.  Examining  Figure  6-3  we  notice  a  deviation  from  the  trend 
at  an  SNR  level  of  5  d8  where  the  standard  deviation  in  SNR  increases 
with  respect  to  the  0  dB  measurement.  This  deviation  may  be  due  to 
numerical  problems  in  the  actual  algorithm  or  in  the  stimulations  soft¬ 
ware.  However,  overall  the  SNR  estimates  appear  quite  reasonable  with  a 
standard  deviation  error  of  less  than  1  dB  for  SNR  levels  greater  than 
or  equal  to  -15  dB. 

6.3  SEMI -STATIC  PERFORMANCE  COMPARISON  BETWEEN  THE  GRADIENT  SEARCH 

ALGORITHMS  AND  OTHER  RELATED  MAP  ALGORITHMS 

In  this  section  we  will  compare  our  chosen  gradient  search  algor¬ 
ithm  against  related  MAP  estimators.  The  comparison  will  be  furtner 
broken  down  into  two  categories.  The  first  will  be  a  comparison  against 


related  gradient  search  algorithms.  The  second  will  be  a  comparison 
against  dynamic  programming  algorithms. 

As  will  be  true  for  the  remaining  simulations  presented  in  this 
chapter,  all  estimates  will  be  from  the  midpoint  of  the  data  block  and 
the  integration  time  T  per  measurement  is  set  to  5  seconds. 
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Figure  6-4  shows  the  RMS  time  delay  difference  errors  versus  SNR 
for  the  following  gradient  search  algorithms. 

1.  The  standard  gradient  search  algorithm  highlighted  in 
Section  6.2.  The  proportional ity  constant  K  =  1.5,  the  block 
size  M  =  10,  the  number  of  lag  points  N  =  9,  and  the  integra¬ 
tion  time  T  =  5  seconds. 

2.  The  standard  gradient  search  ML  estimator  followed  by  a  sepa¬ 
rate  quadratic  LMSF  smoother.  For  this  algorithm  there  is  no 
a  priori  information  penalty  function.  The  block  size  M  =  10, 
the  number  of  lag  points  N  =  9,  and  the  integration  time 

T  =  5  seconds. 

3.  The  standard  gradient  search  algorithm  with  prewhitening. 
(Identical  parameters  as  the  standard  gradient  search  algor¬ 
ithm.  ) 
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The  above  algorithms  are  discussed  in  Chapters  3  3nd  Chapter  1. 


Examining  Figure  6-4  we  can  make  the  following  observations. 

1.  The  gradient  search  algorithm  with  prewhitening  is  at  best 
equivalent  to  the  nominal  gradient  search  algorithm.  Since 
the  prewhitening  involves  a  lot  more  processing  for  no  gain, 
we  can  conclude  that  prewhitening  is  not  a  viable  option. 

2.  The  gradient  ML  estimator  followed  by  the  LMSF  smoother  demon¬ 
strates  instability  at  SNR  levels  below  -18  dB.  Note  that  the 
-20  dB  point  represents  one  successful  estimation  in  five 
attempts  (i.e.,  four  attempts  diverged).  This  experiment 
verifies  our  conclusion  in  Chapter  3.  An  ML  estimation  over  a 
block  of  data  followed  by  an  LMSF  smoother  is  not  equivalent  to 
a  simultaneous  ML  estimation  and  LMSF  a  priori  information 
penalty  function  operating  on  a  block  of  data. 

Figure  6-5  shows  the  RMS  time  delay  difference  error  versus  SNR 
for  the  gradient  search  algorithm  and  the  following  dynamic  programming 
algorithms. 

1.  A  dynamic  programming  algorithm  with  a  Sum  of  Square  Residual 
(SSR)  a  priori  information  penalty  function.  (This  is  the 
same  a  priori  information  penalty  function  used  by  the 
gradient  search  algorithm.) 


2.  A  dynamic  programming  algorithm  with  a  Sum  of  Squared  Differ¬ 
ences  (SSO)  a  priori  information  penalty  function.  (This  is 


the  a  priori  information  penalty  function  derived  by 
Wolcin.6,7,8) 

3.  A  dynamic  programming  algorithm  with  an  SSR  a  priori  informa¬ 
tion  penalty  function  which  contains  a  spectral  estimation 
algorithm  within  the  ML  estimator. 

Note  all  algorithms  used  the  optimum  parameters  found.  Block  size 
M  =  16,  proportional ity  constant  K  =  .075,  number  of  lag  points  N  »  9, 
and  integration  time  T  =  5  seconds. 

Examining  Figure  6-5  we  can  make  the  following  observations. 

1.  Our  first  observation  is  that  the  dynamic  programming  algor¬ 
ithm  containing  the  spectral  estimation  parallels  the  dynamic 
programming  SSR  algorithm  until  -10  dB.  At  SNR  levels  below 
-10  dB  the  spectral  dynamic  programming  algorithm  diverges. 

The  source  of  the  instability  is  apparently  the  spectral 
estimation  which  we  shall  examine  in  more  detail  in  a  follow¬ 
ing  section. 

2.  The  second  observation  is  that  although  the  dynamic  program¬ 
ming  SSD  has  better  performance  at  high  SNR  it  degrades  in  per 
formance  with  respect  to  the  dynamic  programming  SSR  algorithm 
at  low  SNR.  The  SSR  algorithm  appears  to  have  a  lower  thresh¬ 
old  than  the  SSD  algorithm;  however,  this  may  not  actually  be 


single  successful  measurement  in  five  attempts.  However,  the 
superior  performance  at  low  SNR  for  the  SSR  algorithm  over  the 
SSO  algorithm  is  our  reason  for  selecting  the  LMSF  a  priori 
information  penalty  function. 

3.  The  third  observation  is  that  the  gradient  search  algorithm  is 
clearly  superior  in  RMS  time  delay  difference  error  and 
threshold  than  any  of  the  dynamic  programming  algorithms. 
However,  the  current  gradient  search  algorithm  requires  at 
least  four  times  the  computation  time  as  the  current  dynamic 
programming  algorithm. 

4.  The- final  observation  is  the  obvious  quantization  problem  at 
high  SNR  for  the  dynamic  programming  algorithms.  This  point 
is  highlighted  by  the  RMS  time  delay  difference  error  being 
smaller  at  0  dB  than  at  5  dB  as  discussed  in  Chapter  5. 

6.4  GRADIENT  SEARCH  ALGORITHM  DYNAMIC  PERFORMANCE 

In  this  section  we  shall  examine  the  degradation  in  the  gradient 
search  algorithm  with  increased  target  dynamics.  We  have  discussed 
previously  the  difficulty  in  performing  precise  realistic  dynamic  test¬ 
ing.  Therefore  we  claim  only  a  relative  measure  of  the  effects  of 
dynamics  on  the  gradient  search  algorithm.  We  shall  also  present  the 
performance  of  the  gradient  search  algorithms  rate  estimate.  The 
acceleration  estimate  for  the  gradient  search  algorithm  did  not  yield 
meaningful  estimates  and  therefore  shall  not  be  presented. 


The  increased  target  dynamics  represents  a  doubling  and  quadru¬ 
pling  of  the  nominal  target  speed.  Although  the  target  speed  has  been 
doubled  twice,  the  average  time  delay  difference  rate  only  increased 


1.  At  high  SNR  (>-5  dB)  the  increased  target  dynamics  results  in 
a  quite  noticeable  degradation  in  RMS  time  delay  difference 
error.  At  lower  SNR  this  difference  in  RMS  time  delay  differ¬ 
ence  error  is  reduced.  This  result  may  be  due  to  initializa¬ 
tion  transients  which  are  masked  at  lower  SNR. 

2.  At  low  SNR  the  RMS  time  delay  difference  errors  for  the  double 
nominal  target  dynamic  scenario  are  effectively  equivalent  to 
the  corresponding  RMS  time  delay  difference  errors  for  the 
nominal  target  dynamic  scenario.  The  apparent  discrepancy  at 
-22  dB  may  be  attributed  to  marginal  measurement  accuracy  at 
very  low  SNR. 

3.  The  RMS  time  delay  difference  error  for  the  quadruple  nominal 
target  dynamic  scenario  has  a  noticeable  increase  in  error  over 
the  RMS  time  delay  difference  error  for  the  lower  target 
dynamic  scenarios.  Apparently,  the  quadruple  nominal  target 
dynamics  is  sufficient  to  degrade  the  gradient  search  algor¬ 
ithm  at  al 1  SNR  levels . 

Figure  6-7  shows  the  RMS  time  delay  difference  rate  error  versus 
SNR  for  the  same  target  dynamic  scenarios  as  in  Figure  6-6.  The  RMS 
time  delay  difference  rate  errors  will  be  in  units  of  fractional  lag 
point  spacing  per  minute.  The  averaging  time  delay  difference  rate  for 
the  nominal  target  dynamic  scenario,  double  nominal  target  dynamic 
scenario  are  1.8,  3.5,  and  6.8  fractional  lag  point  spacing  per  minute, 
respectively. 
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Examining  Figure  6-7  we  can  make  the  following  observations  about 
the  rate  estimate  of  the  gradient  search  algorithm.  Note,  however,  that 
the  time  delay  difference  rate  estimate  results  from  the  coefficients 
of  the  quadratic  LMSF  a  priori  information  penalty  function  (Equa¬ 
tion  3-106). 

1.  Our  first  observation  is  that  the  maximum  RMS  time  delay 
difference  rate  errors  are  significantly  below  the  average 
time  delay  difference  rate  for  all  three  target  dynamic  sce¬ 
narios  at  all  SNR  levels  measured.  This  result  states  that  the 
time  delay  difference  rate  estimate  is  a  valuable  and  meaning¬ 
ful  estimate  at  all  SNR  levels. 

2.  At  high  SNR  we  notice  that  doubling  and  quadrupling  the  target 
dynamics  (actual  1.95  and  3.8)  has  increased  our  RMS  time 
delay  difference  rate  error  by  factors  of  2.04  and  4.74.  This 
result  indicates  that  the  RMS  time  delay  difference  rate  error 
is  proportional  to  the  actual  time  delay  difference  rate  with 
an  increased  degradation  as  the  actual  time  delay  difference 
rate  becomes  larger. 
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3.  Our  most  striking  observation  is  the  flat  RMS  time  delay 

difference  rate  error  versus  SNR  at  high  SNR  for  all  three  tar- 


6.5  GRADIENT  SEARCH  ALGORITHM  PERFORMANCE  WITH  SPECTRAL  MISMATCH 


The  gradient  search  algorithm  was  derived  under  the  assumption  of 
a  flat  target  power  spectrum.  In  this  section  we  shall  examine  the 
degradation  in  the  gradient  search  algorithm  RMS  time  delay  difference 
error  with  respect  to  target  power  spectrums  not  conforming  to  our 
model.  Our  target  power  spectrums  will  be  limited  to  linear  models 
(i.e.,  the  spectrum  power  changes  linearly  with  frequency  from  the 
lower  cut  off  frequency  f^  to  the  upper  cut  off  frequency  f2).  The 
following  is  a  description  of  the  power  spectral  models  used  in  the  sim¬ 
ulation  results.  Note  that  although  the  power  spectral  models  differ 
they  have  been  normalized  to  contain  equivalent  total  power  over  the 
band  pass  (fj^  through  f2). 

Model  1  10  log  (f2/f x )  3  -6  dB 

Model  2  10  log  (f2/fj)  =  -3  dB 

Model  3  10  log  (f2/f1)  =  0  dB 

Model  4  10  log  (f2/f1)  =  3  dB 

Model  5  10  log  (f 2/f x )  =  6  dB 

Although  the  power  spectral  models  are  specified  in  dB,  they  are  linear 
in  power  and  not  in  dB. 

Figure  6-8  show  the  RMS  time  delay  difference  errors  versus  SNR 
for  the  five  target  power  spectral  models  outlined  above.  The  algorithm 
under  test  is  the  gradient  search  algorithm  highlighted  in  Section  6.2. 


Search  (K 


The  parameters  are:  the  number  of  lag  points  N  =  9;  the  block  size  M  * 
10;  the  proportional ity  constant  K  =  1.5;  and  the  integration  time 
T  =  5  seconds. 

Examining  Figure  6-8  we  can  make  the  following  observations. 

1.  From  an  SNR  level  of  -5  dB  to  -22  dB  the  RMS  time  delay  differ¬ 
ence  error  among  the  various  target  power  spectral  models  do 
not  differ  by  much  more  than  2  dB.  This  result  states  that  our 
gradient  search  estimation  algorithm  is  not  very  sensitive  to 
moderate  target  power  spectral  modeling  mismatches.  Since  in 
practice  target  power  spectral  mismatches  are  quite  likely,  it 
is  fortunate  that  our  estimator  is  robust  with  respect  to 
these  modeling  errors. 

2.  From  about  -10  dB  to  -18  dB  the  RMS  time  delay  difference 
error  for  the  positive  power  spectral  slope  models  (models  4 
and  5)  are  less  than  the  RMS  time  delay  difference  error  for 
the  no  mismatch  zero  power  spectral  slope  (model  3).  This 
result  is  not  surprising  since  the  positive  power  spectral 
slope  emphasizes  the  higher  frequencies.  Substituting  Equa¬ 
tion  3-92  into  Equation  6-1  and  performing  numerical  integra¬ 
tion,  we  can  predict  the  expected  gain  or  loss  in  RMS  time 
delay  difference  error  due  to  our  target  power  spectrum 
models.  The  results  of  this  prediction  indicate  about  a  1  dB 
gain  for  the  positive  sloped  power  spectrum  (model  5)  and  a 

1  dB  loss  for  the  negative  slope  power  spectrum  (model  1)  for 


all  SNR  levels.  This  is  in  agreement  with  our  first  observa¬ 
tion.  Therefore  below  -5  dB  the  differences  in  RMS  time  delay 
differences  error  curves  can  be  related  to  theory  relating  to 
the  power  spectral  shapes  rather  than  modeling  mismatches. 

3.  For  SNR  levels  above  -5  dB  there  is  a  distinct  increase  in  RMS 
time  delay  difference  errors  directly  related  to  the  propor¬ 
tion  of  the  power  spectral  model  mismatch  (i.e.,  the  bigger 
the  power  spectral  mismatch  the  greater  the  RMS  time  delay 
difference  error).  This  fact  indicates  that  power  spectral 
modeling  mismatches  (at  least  for  moderate  modeling  mis¬ 
matches)  are  only  a  factor  at  high  SNR.  At  low  SNR  the  noise 
'field  dominates  most  of  the  distinguishable  target  power  spec¬ 
tral  features,  thus  masking  modeling  errors. 

Figure  6-9  presents  the  RMS  spectral  slope  parameter  y  (Equa¬ 
tion  3-92A)  error  versus  SNR.  The  spectral  slope  parameters  were 
estimated  during  the  same  simulation  runs  shown  in  Figure  6-8.  The 
following  is  a  list  of  the  actual  values  of  the  spectral  slope  param¬ 
eter  y  for  each  of  the  five  target  power  spectral  models. 

Model  1  y  =  -0.598 

Model  2  y  =  -0.332 

Model  3  y  =  0 

Model  4  y  =  0.332 

Y  =  0.598 


Model  5 


Examining  Figure  6-9  we  can  make  the  following  observations. 

1.  The  most  important  observation  is  that  below  -10  dB  (-15  dB 
for  certain)  the  spectral  parameter  estimate  is  unstable  and 
meaningless.  This  is  convincing  evidence  why  in  Chapter  3  we 
did  not  want  to  include  a  spectral  parameter  estimate  as  an 
integral  part  of  our  ML  estimator.  Note  that  the  dynamic  pro¬ 
gramming  ML  estimator  which  contained  a  spectral  parameter 
estimation  as  an  integral  part  diverged  at  SNR  levels  below 
-10  dB  (see  Figure  6-5).  From  Figure  6-8  we  conclude  that 
modeling  mismatches  tend  to  be  masked  by  noise  below  a  -5  dB 
SNR  level.  In  support  of  that  conclusion  we  now  conclude  that 
our  spectral  shape  estimation  is  meaningless  below  -10  dB. 

2.  Above  -10  dB  the  RMS  spectral  shape  estimates  errors  are 
reasonable  proportioned  to  their  actual  values,  thus  we  have 
acceptable  estimators. 

3.  The  larger  power  spectrum  model  mismatch,  the  greater  is  the 
RMS  error  in  the  spectral  shape  parameter.  At  higher  SNR  this 
degradation  in  spectral  shape  parameter  estimation  becomes 
larger.  At  higher  SNR,  target  power  spectrum  modeling  mis¬ 
matches  are  a  more  significant  factor  on  the  algorithm. 
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6.6  INITIALIZATION  PERFORMANCE 

Most  theoretical  work  and  corresponding  empirical  observations 
concern  themselves  with  the,  steady  state  performance  of  the  estimator 
or  algorithm  in  question.  In  Section  6.2  through  6.5  we  concerned  our¬ 
selves  with  the  steady  state  performance  of  the  gradient  search  and 
related  algorithms.  We  have  advertised  that  our  recursive  batch  algor¬ 
ithms  will  have  good  initial ization  properties.  In  this  section  we 
shall  empirically  examine  the  initialization  properties  of  the  dynamic 
programming  (SSR)  and  gradient  search  algorithms. 

Many  algorithms  have  good  initialization  characteristics  for 
unrealistic  zero  target  dynamics.  We  shall  examine  initialization 
characteristics  of  our  algorithms  with  both  target  dynamics  and  initial 
estimate  errors.  Since  initialization  is  a  stochastic  process,  our 
unit  of  measure  shall  be  the  number  of  successful  initializations  for 
given  initialization  conditions  in  10  random  attempts.  We  shall  per¬ 
form  initialization  tests  on  two  target  dynamic  scenarios.  The  first 
will  be  the  nominal  target  dynamic  scenario  illustrated  in  Figure  5-2 
and  the  second  will  be  the  double  nominal  target  dynamic  scenario.  Thi 
scenario  results  from  doubling  the  target  speed  in  the  nominal  target 
dynamic  scenario.  Although  the  average  time  delay  difference  rate  does 
not  exactly  double  (1.95)  due  to  nominal  observer  motion,  certainly  at 
the  critical  initialization  stage  the  time  delay  difference  rate  should 
be  almost  exactly  doubled. 
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Using  both  dynamic  target  scenarios,  we  shall  initialize  our 
algorithms  under  test  -2,  -1,  0,  1,  and  2  lag  points  away  from  the 
actual  initial  time  delay  difference.  We  shall  perform  this  experiment 
for  three  SNR  levels  (-15  dB,  -18  dB  and  -20  dB). 

The  first  algorithm  that  we  shall  examine  is  the  dynamic  program¬ 
ming  (SSR)  algorithm  described  in  detail  in  Chapter  4.  Since  this 
algorithm  simultaneously  examines  a  block  of  potential  time  delay 
difference  trajectories,  we  have  claimed  that  it  should  have  very  good 
initialization  properties. 

The  following  is  a  list  of  the  choosen  parameters  for  the  dynamic 
programming  (SSR)  algorithm. 

Block  size  M  =  16 

Number  of  lag  points  N  =  9 

Proportionality  constant  K  =  .075 

Integration  time  T  =  5  seconds 

Number  of  subdivisions  per  lag  spacing  NP  =  4 

Maximum  subdivision  change  per  measurement  update  NPMAX  =  8. 

The  results  for  the  initialization  simulations  for  the  dynamic 
programming  SSR  algorithm  are  given  in  Table  6-1  for  the  nominal  target 
dynamic  scenario  and  in  Table  6-2  for  the  double  nominal  target  dynamic 
scenario. 


Table  6-1.  Dynamic  Programming  (SSR)  Nominal  Target  Dynamics 


Initial 

Laq  Point 

Offset 

SNR 

-2 

-1 

0 

+1 

+2 

-15 

10/10 

10/10 

10/10 

10/10 

10/10 

-18 

9/10 

9/10 

9/10 

10/10 

10/10 

-20 

2/10 

4/10 

2/10 

4/10 

6/10 

Ratio  of  successful  initializations  to 
attempts 


Table  6-2.  Dynamic  Programming  (SSR)  Double 
Nominal  Target  Dynamics 


Initial  Lag  Point  Offset 


SNR 

(dBJ 

-2 

-1 

0 

+1 

+2 

-15 

10/10 

10/10 

10/10 

10/10 

10/10 

-18 

10/10 

9/10 

9/10 

9/10 

10/10 

-20 

0/10 

3/10 

4/10 

5/10 

3/10 

Ratio  of  successful  initializations  to 
attempts 


Examining  Tables  6-1  and  6-2  we  can  make  the  following  observa¬ 
tions  on  the  dynamic  programming  SSR  algorithm  initialization  perform¬ 
ance. 

1.  Our  first  observation  is  that  the  dynamic  programming  (SSR) 
algorithm  has  excellent  initialization  properties  down  to 
-18  dB.  The  increased  target  dynamics  appears  to  have  little 
affect  at  SNR  levels  of  -15  dB  and  -18  dB. 

2.  At  -20  dB  the  apparent  initialization  performance  falls  con¬ 
siderably  off.  However,  -20  dB  is  near  the  steady  state 
threshold  for  the  dynamic  programming  (SSR)  algorithm.  Many 
of  the  failed  initialization  attempts  occurred  at  times'mary 
minutes  after  the  initialization  and  potentially  could  repre¬ 
sent  unfortunate  noise  field  variations. 

3.  As  would  be  expected,  at  -20  dB  the  double  nominal  target 
dynamics  scenarios  have  poorer  initialization  ratios  than  the 
nominal  target  dynamic  scenario. 

4.  The  initialization  success  ratios  for  the  initial  positive  lag 
point  offset  are  in  general  superior  to  the  0  or  negative  ini¬ 
tial  lag  point  offsets.  Since  both  dynamic  target  scenarios 
are  progressing  in  a  positive  lag  point  sense,  this  fact  indi¬ 
cates  that  better  performance  can  be  obtained  by  initializing 
at  a  time  delay  difference  point  leading  the  dynamic  target. 
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Overall,  we  claim  that  the  dynamic  programming  (SSR)  algorithm  has 
excellent  initialization  properties. 

We  shall  now  examine  the  gradient  search  algorithm  described  in 
Chapter  4  and  examined  in  Sections  6.1  through  6.5.  We  have  tried  to 
give  this  algorithm  as  good  initialization  properties  as  the  dynamic 
programming  SSR  algorithm.  However,  the  gradient  search  algorithm  is 
more  dependant  on  the  initial  estimate  and  does  not  perform  multiple 
simultaneous  searches. 

The  parameter  choices  for  the  gradient  search  algorithm  are 
repeated  below  for  convenience. 

Block  size  M  =  10 

Number  of  lag  points  N  =  9 

★ 

Proportionality  constant  K  *  1.5 

Averaging  time  T  =  5  seconds. 

Simulations  have  indicated  that  although  larger  values  of  the  pro¬ 
portionality  constant  K  yield  better  steady  state  estimator  perform¬ 
ance,  the  initialization  performance  degrades  with  higher  proportion¬ 
ality  constant  values.  This  point  is  illustrated  in  Figures  6-10 
through  6-12.  They  each  plot  the  time  delay  difference  error  versus 
time  for  the  gradient  search  algorithm  initialized  on  target  (0  lag 

*  The  final  value  of  the  proportional ity  constant  K  will  be  1.5. 


offset)  for  a  SNR  level  of  -18  dB.  However,  in  Figure  6-10  the  propor¬ 
tionality  constant  is  set  to  1.5;  in  Figure  6-11  the  proportionality 
constant  K  is  set  to  0.5;  and  in  Figure  6-12  the  proportionality  con¬ 
stant  is  set  to  0.1. 

Examining  Figure  6-10  (K  =  1.5),  we  can  observe  a  large  initiali¬ 
zation  transient.  In  several  other  simulations  under  similar  condi¬ 
tions,  this  initialization  transient  has  led  to  a  loss  of  track  or 
instabi 1 ity. 


Examining  Figures  6-11  and  6-12  for  the  lower  value  of  the  propor¬ 
tionality  constant  K,  we  can  see  that  the  initialization  transient  is 
removed  but  the  random  variations  are  greater.  Similar  simulations  at 
-20  dB  have  shown  poor  transient  behavior  for  the  K  *  0.5  proportion¬ 
ality  constant  with  only  K  s  0.1  proportional ity  constant  having  satis¬ 
factory  characteristics. 


The  basic  problem  is  that  increasing  the  proportional ity  constant 
increases  the  weight  of  the  LMSF  a  priori  information  function.  Our 
a  priori  information  is  that  the  target  time  delay  difference  traject¬ 
ory  is  a  quadratic  function  in  time  over  the  block  size.  The  a  priori 
information  depends  largely  on  the  actual  ML  estimates.  In  the  case  of 
initialization,  the  initial  estimate  and  the  following  transient  ML 
estimates  are  unreliable  for  use  by  the  LMSF  a  priori  information  func¬ 
tion.  Therefore,  relying  too  much  initially  on  the  LMSF  a  priori 
information  function  compounds  the  error  which  is  difficult  to  correct 
for  as  time  goes  on. 
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We  do,  however,  desire  to  have  the  optimum  smoothing  qualities  of 
the  larger  proportionality  constant  K  once  in  steady  state.  To  obtain 
the  best  of  both  initialization  and  steady  state  properties,  the 
following  initialization  algorithm  was  adopted. 

1.  From  time  =  0  seconds  until  time  =  110  seconds  (12  updates 
after  the  block  is  filled),  the  proportionality  constant  is  set 
to  0.1. 

2.  From  time  *  115  seconds  until  time  =  250  seconds  (28  updates), 
the  proportionality  constant  is  incremental ity  update  by  0.05. 

3.  From  time  255  seconds  onwards,  the  proportional  ity  constant  !< 
is  fixed  at  1.5. 

The  slow  transition  in  step  2  from  K  =  0.1  to  K  =  1.5  is  required 
to  minimize  introducing  transients  into  the  algorithm. 

Figure  6-13,  performed  on  the  identical  scenario  as  Figures  6-10, 
6-11,  and  6-12,  illustrates  the  performance  of  our  new  initialization 
algorithm. 

Examining  Figure  6-13,  we  can  see  at  the  beginning  the  same  output 
as  Figure  6-12  for  the  K  =  0.1  proportional ity  constant.  At  the  end. 
Figure  6-13  resembles  Figure  5-10  for  the  K  =  1.5  proportional ity  con¬ 
stant.  In  between,  there  are  no  apparent  initial  izaiton  transients. 


We  do  not  claim  that  our  initialization  algorithm  is  optimum.  It 
is  certainly  a  subject  for  future  work.  However,  it  is  the  algorithm 
for  which  we  present  the  equivalent  initialization  test  results  for  the 
gradient  search  algorithm  as  we  did  earlier  for  the  dynamic  programming 
(SSR)  algorithm. 


Tables  6-3  presents  the  nominal  target  dynamic  initial izaiton  test 
results  for  the  gradient  search  algorithm  and  Table  6-4  presents  the 
double  nominal  target  dynamics  initialization  test  results  for  the 
gradient  search  algorithm. 


Examining  Tables  6-3  and  6-4  we  can  make  the  following  observa¬ 
tions  about  the  initialization  properties  of  the  gradient  search  algor¬ 
ithm. 


1.  At  -15  dB  for  either  target  dynamic  scenario,  the  gradient 
search  algorithm  is  perfect  in  all  initialization  attempts. 


2.  At  -18  dB  we  notice  a  slight  degradation  in  successful  ini¬ 
tializations.  However,  at  double  target  dynamics  this  degra¬ 
dation  becomes  significant  at  an  initial  offset  of  +2  lag 
points. 


3.  At  -20  dB  for  the  nominal  target  dynamics,  the  initialization 
performance  is  excellent  for  negative  and  zero  initial  lag 
point  offsets  but  poor  at  positive  initial  lag  point  offsets. 
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4.  At  -20  dB  for  the  double  normal  target  dynamics,  the  initial¬ 
ization  falls  way  off.  Once  again,  the  negative  and  zero  ini¬ 
tial  lag  point  offsets  have  better  initialization  performance 
than  the  positive  lag  point  offsets. 

5.  It  is  apparent  that  the  gradient  search  algorithm  initializes 
significantly  better  for  a  lagging  or  negative  initial  lag 
poirt  offsets  than  for  leading  or  positive  initial  lag  point 
offsets.  This  result  is  the  opposite  of  the  dynamic  program¬ 
ming  (SSR)  algorithm.  The  reason  for  the  above  is  almost 
certainly  the  result  of  the  gradient  search  algorithm  ini¬ 
tially  estimating  the  target  time  delay  difference  trajectory 
slope  in  the  wrong  direction.  (A  lead  offset  will  result  in  a 
negative  initial  time  delay  difference  slope  when  the  actual 
time  delay  difference  trajectory  slope  is  positive.)  Although 
steps  have  been  taken  to  minimize  this  problem,  they  obviously 
have  not  been  totally  successful.  Therefore,  there  is  obvious 
improvements  which  can  be  incorporated  in  the  gradient  search 
algorithm  initialization  procedure. 

Overall,  the  gradient  search  initialization  characteristics  are 
slightly  worse  than  the  dynamic  programming  (SSR)  algorithm  initializa¬ 
tion  characteristics.  However,  the  gradient  search  algorithm  is  capabl 
of  operating  at  SNR  level  at  least  3  dB  lower  than  the  dynamic  program¬ 
ming  (SSR)  algorithm.  Therefore,  we  can  conclude  that  the  gradient 
search  initial izaiton  procedure  is  not  as  good  as  the  dynamic  program¬ 
ming  initialization  procedure  and  that  there  is  a  definite  gap  between 


CHAPTER  7 


SUMMARY  CONCLUSIONS,  AND  RECOMMENDATIONS 

We  have  studied  in  detail  the  problem  of  estimating  time  delay 
difference  parameters  of  a  moving  broad  spectrum  target  from  two  obser¬ 
vation  sensors.  The  objective  of  this  study  was  to  derive  recursive 
batch  estimation  algorithms  with  good  initialization,  threshold  and 
dynamic  properties.  The  practical  problem  which  motivated  our  study  is 
primarily  in,  but  not  limited  to,  the  field  of  underwater  acoustics. 

The  acoustical  problem  is  to  estimate  the  time  delay  difference  and/or 
bearing  trajectory  of  a  moving  broad  spectrum  acoustical  target. 

Our  theoretical  approach  to  the  above  problem  was  to  derive  a  MAP 
estimator  similar  to  the  estimator  derived  by  Wolcin6  for  narrowband 
signal  parameters  estimation.  The  MAP  approach  was  choosen  to  take 
advantage  of  known  a  priori  information  on  expected  acoustical  target 
trajectories.  We  assumed  that  over  a  significant  period  of  time,  time 
delay  difference  trajectories  will  be  smooth,  continuous  and  can  be 
modeled  by  a  low  order  (quadratic)  polynomial  in  time. 

The  derivation  of  the  maximum  likelihood  (ML)  portion  of  the  MAP 

2  5 

estimator  was  along  similar  lines  adopted  by  Knapp  and  Carter  and  Ng. 
Our  principal  contribution  in  this  area  was  the  precised  derivation  of 
the  expected  cross  correlation  (MISMATCH)  function  in  the  presence  of  a 
linear  time  delay  difference  trajectory.  This  derivation  contained  an 
additional  factor  predicting  a  coherency  loss  with  increase  observation 
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time  and  time  delay  difference  rate.  In  addition,  we  have  expanded  the 
derivation  of  the  mismatch  function  for  the  more  general  case  of  a 
linear  target  power  spectrum.  This  work  includes  an  estimator  for  a 
parameter  related  to  the  slope  of  the  linear  target  power  spectrum. 

Our  derived  ML  estimator  assumes  as  input  the  discrete  noisy 
measurements  of  the  cross  correlation  between  our  two  observation  sens¬ 
ors  averaged  for  T  seconds  and  evaluated  at  N  equally  spaced  lag  delays 
centered  at  zero  lag  delay.  The  noise  components  between  measured  lag 
points  are  assumed  independent,  although  the  correlated  noise  compon¬ 
ents  case  is  examined. 

A  block  of  input  data  consists  of  a  fixed  number  of  M  consecutive 
measured  lag  point  data  updates.  A  recursive  batch  process  is  a  batch 
process  where  at  every  update  the  oldest  measured  data  is  deleted  from 
the  data  block  and  the  newest  measured  data  is  added  to  the  data  block. 
After  every  update,  the  recursive  batch  process  calculates  new  esti¬ 
mates  based  solely  on  the  new  block  of  data. 

5 

Similar  to  Ng,  we  adopted  two  different  time  delay  difference 
trajectory  models  to  the  same  problem  to  facilitate  our  derivation  of  a 
batch  estimator.  To  maximize  response  time,  minimize  coherency  loss 
and  provide  a  tractable  model  for  our  ML  estimator,  we  have  assumed  that 
the  unknown  time  delay  difference  trajectory  is  a  piecewise  linear 
function  over  the  block  of  measurements  of  interest.  The  piecewise 
inear  segments  correspond  to  individual  measurement  update  or  observa¬ 
tion  times.  In  order  to  gain  the  advantage  of  long  observation  time,  we 
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have  assumed  that  our  unknown  time  delay  difference  trajectory  can  be 
modeled  by  a  quadratic  function  in  time  over  the  block  of  time  of  inter 
est.  To  tie  these  two  models  together,  we  have  derived  an  a  priori 
information  penalty  function  based  on  minimizing  the  squared  error  of 
the  individual  ML  time  delay  difference  estimate  over  our  block  of  data 
measurements  with  respect  to  a  quadratic  time  delay  difference  model. 
The  concept  of  this  LMSF  a  priori  information  penalty  function  operat¬ 
ing  simultaneously  with  our  individual  ML  estimates  is  unique  from 
either  the  work  of  Ng^  or  Wolcin.*’ 

From  our  theoretical  work  we  have  derived  two  major  algorithmic 
implementations.  The  first  is  based  on  a  dynamic  programming  technique 
adopted  by  Wolcin.7  We  claim  credit  for  obtaining  an  efficient  recur¬ 
sive  implementation  of  this  algorithm.  In  our  recursive  implementation 
of  the  LMSF  a  priori  information  function,  we  have  used  a  technique 
developed  by  Ng.^  However  we  have  uncovered  potential  numerical 
instability  in  Ng's  algorithm  and  have  derived  a  numerical  technique  to 
eliminate  this  instability. 

The  second  algorithmic  implementation  is  based  on  numerical  grad¬ 
ient  search  techniques.  Although  the  numerical  techniques  are 
straightforward,  neither  Wolcin,  Ng,  nor  any  known  researcher  has 
applied  these  recursive  search  techniques  to  this  type  of  problem. 


Other  theoretically  related  algorithms  were  also  implemented  for 
comparison  studies.  The  most  notable  of  these  is  a  algorithm  which 
prewhitens  the  correlated  input  correligram  lag  point  data  to  effec¬ 
tively  obtain  independent  noise  on  all  the  input  correlogram  lag 
points. 

Using  empirical  result  based  on  simulations  and  where  possible 
theoretical  predictions,  we  have  sought  to  find  the  optimum  algorithmic 
parameter  selections  for  our  two  major  time  delay  estimation  algor¬ 
ithms.  The  algorithmic  parameters  we  have  researched  are  listed  below. 

1.  The  proportionality  constant  K  between  the  ML  estimator  and 
the  a  priori  information  penalty  function. 

2.  The  number  M  of  data  measurements  in  our  data  block  and  the 
integration  time  T  per  measurement. 

3.  The  number  of  lag  points  N  calculated  for  each  measurement 
upd ate. 

4.  The  lag  spacing  Tau  between  measured  lag  points. 

The  following  major  conclusions  were  drawn  from  the  above  simula¬ 
tion  experiments.  The  simulations  were  for  the  most  part  performed 
using  a  moderate  dynamic  time  delay  difference  trajectory.  The  stimu¬ 
lator  used  for  the  simulation  experiments  was  independently  developed. 
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The  stimulator  produced  simulated  correlogram  lag  point  data  with 
realistic  correlated  noise  added  to  the  deterministic  lag  points. 


1.  There  are  optimal  values  for  the  proportional ity  constant  K 
for  both  the  gradient  search  and  dynamic  programming  algor¬ 
ithms.  The  optimum  choice  of  K  is  a  function  of  SNR,  block 
size  M,  and  target  dynamics.  A  corollary  of  this  observation 
is  that  for  reasonable  target  motion  the  proper  introduction 
of  a  priori  information  (K  >  0)  will  always  improve  the  vari¬ 
ance  of  our  estimators. 


2.  For  the  dynamic  programming  algorithm,  the  ML  estimator 

(K  =  0)  has  little,  if  any,  dependence  on  block  size  after  a 
minimum  block  size  of  4.  The  minimum  dependence  of  block  sizes 
of  4  or  less  is  due  to  the  weak  piecewise  linear  assumption  of 
the  ML  estimator.  For  the  moderate  test  dynamics,  the  ML 
estimator  is  always  improved  by  moderate  increased  averaging 
time. 


The  above  observation  should  also  be  true  for  the  gradient 
search  algorithm  as  well.  However,  numerical  and  algorithmic 
problems,  especially  at  high  SNR,  indicate  an  optimum  choice 
of  block  size. 


3.  With  the  introduction  of  optimum  a  priori  information,  there  is 
an  optimum  choice  of  block  size  dependent  on  integration  time 
and  target  dynamics.  Barring  numerical  problems  due  to  the 
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increased  number  of  unknowns,  the  following  is  the  explanation 
for  the  dependence  of  block  size  on  averaging  time  and  target 
dynamics.  The  more  numbers  of  measurements  and  the  more 
severe  the  target  dynamics  the  greater  the  probability  of  vio¬ 
lating  algorithmic  assumptions  (quadratic  time  delay  differ¬ 
ence  trajectory  model ).  The  longer  the  integration  time  per 
measurement  the  more  likely  a  modeling  mismatch  will  occur. 
(For  large  block  sizes,  shorter  integration  times  actually  can 
improve  algorithm  performance.) 


4.  For  relatively  low  dynamic  targets  the  number  of  lag  points  N 
does  not  affect  algorithmic  performance.  This  conclusion  is 
based  upon  the  realistic  assumption  of  correlated  noise 
between  lag  point  measurements.  If  the  noise  was  independent 
this  conclusion  would  be  valid  after  9  lag  points. 


5.  The  prewhitening  option  does  not  significantly  improve  our  ML 
estimator  in  the  case  of  correlated  noise. 


6.  For  high  target  dynamics  a  small  number  of  lag  points  (3)  may 
degrade  performance.  Since  increasing  the  number  of  lag 
points  increases  computations  for  our  algorithm,  there  exists 
an  optimum  choice  for  the  number  of  lag  points.  This  optimum 
choice  lies  in  the  range  of  3  to  9  lag  point  measurements. 


7.  Given  target  dynamics,  there  is  an  optimum  choice  of  averaging 
time  T  for  the  measured  lag  point  data  supplied  to  our  ML 


estimator.  The  more  severe  the  target  time  delay  rate,  the 
shorter  the  optimum  averaging  time  T.  Theoretically,  this 
optimal  is  a  balance  between  averaging  time  and  coherency  loss 
due  to  target  motion. 

8.  Theoretically,  the  choice  of  lag  point  spacing  will  not  affect 
our  ML  estimators  performance  unless  it  is  larger  than  the 
Nyquist  rate  for  the  highest  processed  frequency. 

In  practice,  too  small  choices  of  lag  point  spacing  severely 
degrade  algorithmic  performance.  Therefore,  the  optimum 
choice  of  lag  point  spacing  is  in  the  order  of  1/2  the  Nyquist 
rate  to  the  actual  Nyquist  rate' for  the  highest  processed 
frequency. 

Our  final  effort  was  to  illustrate  the  performance  of  the  chosen 
gradient  search  algorithm  with  respect  to  a  basic  standard,  other 
related  algorithm  alternatives,  and  varying  spectral  and  dynamic  test 
conditions.  The  basic  reference  was  selected  as  the  Cramer-Rao  Lower 
Bound  (CRLB)  for  coherent  integration  of  a  stationary  target  (no 
a  priori  information).  We  shall  name  this  basic  reference  the  Maximum 
Likelihood  Estimate  (MLE)  CRLB.  The  performance  tests  consisted  of 
Monte  Carlo  simulations  conducted  on  the  gradient  search  and  other 
related  algorithms.  Once  again,  our  stimulator  assumes  correlated 
noise  between  lag  point  measurements. 


estimation  performance  could  be  attributed  to  the  actual  tar¬ 
get  power  spectral  shape  (i.e.,  not  to  our  incorrect  modeling 
of  the  target  power  spectral  shape). 

The  separate  estimation  of  a  spectral  parameter  was  only 
stable  for  SNR  levels  of  -10  d8  or  greater.  At  these  SNR 
levels  the  spectral  parameter  estimates  were  reasonable. 

7.  The  initialization  properties  of  the  dynamic  programming 
algorithm  are  excellent  for  moderate  target  dynamics.  The 
dynamic  programming  algorithm  had  good  initialization  proper¬ 
ties  to  within  a  few  dB  of  its  steady  state  performance  limit. 

8.  The  initialization  levels  of  the  gradient  search  algorithm  are 
almost  as  good  as  the  dynamic  programming  algorithm.  However, 
its  steady  state  threshold  is  3  dB  lower  than  the  dynamic 
programming  algorithm.  Therefore,  there  is  a  significant  gap 
between  initialization  threshold  and  steady  state  threshold. 
This  difference  is  increased  in  the  presence  of  increased  tar¬ 
get  dynamics. 

7.1  FUTURE  WORK 

Our  work  has  been  a  complete  research  on  the  development  of 
efficient  recursive  batch  time  delay  difference  estimation  algorithm. 
Two  theoretical  areas  requiring  addition  work  are:  1)  A  theoretical 
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development  on  the  functional  dependencies  of  the  proportionality  con¬ 
stant  K;  and  2)  The  related  theoretical  development  of  a  lower  perfor¬ 
mance  bound  for  the  ML  estimator  and  a  priori  information. 
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The  following  areas  are  recommended  for  future  work  in  the  algor¬ 
ithmic  area. 
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1.  Development  of  an  improved  initialization  scheme  and  related 
theory. 

2.  Development  of  adaptive  parameter  selection.  Optimum  param¬ 
eter  selection  is  related  to  the  target  dynamic  state  and  SNR. 
Therefore,  an  adaptive  parameter  algorithm  could  be  developed 
to  conform  to  changing  target  dynamic  conditions  and  SNR 
levels. 

3.  Extension  of  algorithms  to  include  the  mu  1 ti -parameter  case 
such  as  multiple  frequency  bands,  multiple  beams,  and  multiple 
target 1 s. 
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APPENDIX  A 

ALTERNATIVE  DERIVATION  OF  THE  MISMATCH  FUNCTION 


The  following  is  a  slight  modification  of  derivation  suggested  by 
Ng13  for  the  evaluation  of  Equation  3-54.  Equation  3-54  is  repeated 
below  for  convenience. 


G$(u>)  Lim(Sinc  [(Sjoj  -  co^] 


jou(Adn+A8{T/2) ) 
e  u  )dw  . 


Sine  [($2  <d 


wk)2J 


( A— 1 ) 


jui(Adn+AB(T/2) 

Let  6  *  MAX  (3j,  S2)  and  assume  that  Gs (o>)e  changes  very 

little  over  the  interval  [ (a>k  -  2tt/T)/B,  (c^  +  2tt/T)/0]  ,  then 
Equation  A-l  can  be  rewritten  as  the  following: 


ju>./B(Adn+A8(T/2))  l 

Gs(w)e  /  Lim(Sinc  [  (Sico  -  oj.  )i]  • 

J  0  *** 

Sine  [(S2oj  -  uk)^])|^  .  (A-2) 

As  t  becomes  large  (i.e.,  T  becomes  large),  the  width  of  the 
Sine ( - )  functions  are  much  smaller  than  B,  the  signal  bandwidth. 
Therefore,  we  can  approximate  the  following  expression  from 
Equation  A-2. 
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From  this  point  onwards  we  shall  be  concerned  with  evaluating  the 
integral  portion  of  expression  A-3. 


The  key  idea  we  shall  exploit  is  Parseval's  theorem  stated  below. 
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where 

f^tJ-M-F^co) 

f2(t)++F2(ai) 

are  Fourier  transform  pairs. 

Note  that  if  we  define  f(t)  as  the  following: 


(A-4) 


f(t)  =  f. 


1 1 1  <  2, 


(A-5) 


0,  | t |  >  j 

then  the  corresponding  Fourier  transform  F(uj)  is  defined  as  the 
fol lowing: 

F(u)  =  S  i  n  c  ( u>T  /  2 ) 
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Furthermore,  we  have  the  following  relationships  for  Fourier 


transforms . 


f(at)  ~  7 If  F(I}  =  TaT  Sinc(^} 


j«ht 
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In  order  to  evaluate  expression  A-4,  let 


F^(-u>)  =  Sine  [8^(co  - 


F2(u)  =  Sine  [02(o>  -  ^)|] 


With  Equations  A-7  and  A-8  in  mind,  the  inverse  Fourier 
of  Equations  A-9  and  A-10  is  given  below. 


F~1{F1(oj)}  =  - - 


f{~ )  e  k  1 
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f(f-)  e 

-1 

F  Votui)}  -  - ; - 


(A-7) 

(A-8) 

(A-9) 

(A-10) 

transform 

(A- 1 1 ) 

( A- 1 2 ) 


v<  v  -j*  ->  n>  '>  jr-v^  jr^  \r*v 


where  f(*)  is  defined  by  Equation  A-5  and  diagrammed  below 


Therefore,  using  Equations  A-ll,  A-12n  and  A-4  we  can  write  the 


following. 


oo 

J~  Sinc[(S1oj  -  <ok)^]  Sinc[(82u  -  = 


ts  i  /  c 

"f&T/Z  sla2r 


l  jcok  (1/S2-1/S1)t 


where 


8  =  min  (8^  82) 


The  left-hand  side  of  Equation  A - 1 3  evaluates  to  the  following 


after  integration  and  algebra. 


(A— 14 ) 


jf  Sinc[u)k  A8T/26] 

where 

S  =  MAX  [$1 ,  b2] 

A8  =  ^2  “ 

Combining  the  result  for  expression  A-14  with  expression  A-2  and 
noting  that  Gs(co)  is  equal  to  zero  for  |cu|  >  B,  we  have  the  following 
expression  for  the  integral  of  expression  A-l. 

ju>k/S(Adg+A8T/2) 

Gs(<,'k/8)  - - p -  Sinc[wk  A8T/2  8]  ,  -SB  <  ^  <  BB  . 

(A-15) 

This  is  the  result  identical  to  Equation  3-64. 


APPENDIX  B 


EVALUATION  OF  MISMATCH  FUNCTION  INTEGRALS 

This  appendix  shows  the  analytical  approximations  used  to  evaluate 
integrals  related  to  the  mismatch  function. 

From  Chapter  3  we  have  the  integral  of  Equation  3-68  to  evaluate 
which  is  repeated  below  for  convenience. 


wc+Bw/2 

f  Sinc(aj|^T)  cos (u>[t  +  ( AdQ  +|&T)])du 

*'u)c-Bw/2 


(B-l) 


The  above  integral  can  be  symbolically  written  as  the  following. 


jl  r 

Bw  I 


a>c+Bw/2 


Sine  [A  u]  cos  [8  aj]du 


(B-2) 


ujc-Bw/2 


where 


a  -  AST 
A  -  — 


(B-3) 


and 


B  =  (t.  +  Ad0)  +  ^ 


(B-4) 


we  shal 1  assume  that 

A  oj  <<  1,  jj  -  Bw/2  <  u)  <  ui  +  Bw/2 


;b-b  ; 
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Applying  the  limits  to  expression  B-9  and  using  algebra  and  trig¬ 
onometry  yields ‘the  following  result. 


- - - r  [B  Bw  cos(B  <j)  )  cos(B  Bw/2) 

3Bw  BJ  c 

-  2B  w c  sin(B  ojc)  sin(B  Bw/2) 

+  (B2(u>2  +  (Bw/2)2)  -  2)  cos(B  coc)  sin(B  Bw/2) 

+  id  Bw  B2  sin(B  oj  )  cos(B  Bw/2)]  .  (8-10) 

Note  that  if  the  target  radial  motion  is  zero,  then  AB  *  0  which 
implies  that  A  =  0  which  implies  that  expression  B-10  is  zero. 

To  summarize  our  result,  the  integral  of  expression  B-l  evaluates 
to  the  sum  of  expression  B-7  and  B-10. 

Also  from  Chapter  3  we  have  to  evaluate  the  integrals  from 
Equation  3-94.  The  first  and  third  integrals  from  Equation  3-94  are 
identical  to  the  integral  of  Equation  B-l  which  we  have  already 
evaluated.  The  second  integral  of  Equation  3-94  is  the  result  of  a 
linear  sloped  power  spectrum.  This  integral  is  written  below  for  con- 


•  • » ■ « • 

. .  4 

Vav 

.a 

■ 


venience. 


cos(u>[t.  +  Adg  +  ^|^])  sin(u^|^)da)  . 


(B-ll) 


u> +Bw/2 


r 


Bw  ABT 

The  above  integral  can  be  symbolically  written  as  the  following 


o)c“Bw/2 


a) +Bw/2 


-r 


cos  [A  01]  sin(B  u)]daj 


Bw  ABT 

where  A  and  B  are  defined  by  Equations  B-3  and  B-4,  respectively. 


(B-12) 


d}c-8w/2 


The  integral  shown  in  expression  B-12  can  be  evaluated  exactly. 
From  a  table  of  integrals  we  obtain  the  following  expression 


2  I  cos  [B  -  A 1  co  COsfB  +  Aim 

Bw  ABT  L"  2[B  -  A]  '  2[B  +  A] 


a>c  +  Bw/2 
u>c  -  Bw/2 


(B-13) 


Applying  the  limits  and  using  algebra  and  trigonometry,  expres¬ 
sion  B-13  evaluates  to  the  following. 


■gf  (Sln(u»c[A  +  Bj )  Sinc(|^[A  +  B])  -  sin(u)c[B  -  A])  Sinc(|*[B  -  A])] 

(8-24) 

Expression  B-14  is  the  analytical  result  of  the  integration  of 
expression  B-12. 


To  completely  specify  expression  B-14,  we  should  give  the  limits 
when  ABT  -*•  0,  t^  +  AdQ  0  (or  equivalently  B  -  A  -*■  0)  and 
ABT  +  t.  +  AdQ  +  0  (or  equivalently  B  +  A  -*■  0). 


Lim  [Expression  B-14]  =  u  [cos(<d.  e)  Sinc^  e)  + 
AgT+0  c  c  6 


Adg+t^O 


(cos(|^aug)  -  Sinc(|^  e)]  Sinc(ojc  e)] 


where  e  =  Adg  +  ti  =  B  -  A 


Lim  [Expression  B-14]  = 
A0T+O 


Adg+t^O 


Lim 

ABT^O 

AdQ+iyO 

ABT+Adg+t  ^-*-0 


[Expression  B-14]  =  - 


sin[u)c  e]  sinc[|^  e]/ABT 


APPENDIX  C 

RECURSIVE  LMSF  ERROR  ACCUMULATION  EXPERIMENT 

This  appendix  describes  the  simulation  experiment  to  verify  the 
accumulation  of  errors  in  the  recursive  LMSF  procedure  as  presented  by 
Ng*®  and  outlined  in  Section  4.2.  In  addition,  the  simulation  experi¬ 
ment  will  demonstrate  that  the  modifications  recommended  by 
18 

LaTourette  effectively  eliminate  the  accumulation  of  error  by  the 
recursive  procedure. 

The  simulation  experiment  was  designed  as  follows: 

1.  The  input  signal  time  series  is  described  by  the  following 
equation. 

Yi  =  A  sin(|2L  t.)  +  r.(0,  var),  ti  =  0,  1,  2,  3  ....  N  (C-l) 
where 

Y.  =  the  input  time  series 
A  *  a  constant  amplitude 
T  =  the  period 

r^'s  =  independent  Gaussian  random  variables  with 
mean  =  0  and  variance  =  var 
N  =  the  number  of  points  in  the  time  series. 


C-l 


The  theoretical  matching  function  is  a  quadratic  equation 
listed  below. 


c2  ti  +  C1  ti  +  c0  (C’2) 

where  C2.  Cp  Cq  are  unknown  coefficients. 

The  time  series  Yi  (Equation  C-l)  is  recursively  applied  to 
Equations  4-13,  4-14,  4-15,  and  4-17  which  are  rewritten  below 
in  a  convenient  format. 


M1)  =  li-ld)  +  Yi  -  Yi-M-1 
?i (2)  *  li.1(2)  +  (M  +  1)  Y .  -  7i(l) 

1  i ( 3 )  »  1^(3)  +  (M  +  l)2  Y1  +  21  i (2)  -  11(1) 


(C-3) 

(C-4) 

(C-5) 

(C-6) 


where 

1 . ( k )  =  the  recursive  coefficient  parameters  k  =  1, 
2,  3 


M  +  1  =  the  fixed  number  of  data  point  in  the  block 
2 

SYi  =  the  recursive  sum  of  squared  input  Y .  within  the 
sliding  window. 


•*.  *,•  ^*  *.*,*.  -".  A  ,*  v*.  .•  .*  /  *. 


The  sum  of  the  squared  residual  SSR  for  the  ith  time  frame  is 
obtained  via  Equation  4-16  and  rewritten  below  in  a  convenient 
format. 

SSR.  =  SY*  -  lT  (HT  H)"1  Lj  (C-7) 

where  the  vector  is  defined  below 

l{  «  [li(l),  1  i  (2) ,  1 i (3)1  (C-8) 

and  the  matrix  H  is  defined  in  Equation  4-9. 

Therefore,  the  recursive  parameters  of  interest  are  1^(1), 
1^2)  and  1^(3)  defined  by  Equations  C-3,  C-4,  C-5  and  the 
recursive  estimate  of  the  sum  of  squared  errors  SSRi  defined 
by  Equation  C-8. 

To  determine  the  accumulated  error  in  our  recursive  parameters 
of  interest,  we  shall  calculate  their  values  using  the  follow¬ 
ing  standard  non-recursive  equations. 

From  Equation  4-12,  1^(1),  1^(2)  and  1^(3)  can  be  calculated 
from  the  following  equation. 


C-4 


where  the  vector  AD  is  defined  as 


J  , 


&  S  [Yi*  Vi’  Vm^ 


(C-10 


Therefore, 


'i'1’  -IXh* 


(C-ll 


k=0 


M 

Tf(2)  *£>  -  U 

k=0 


(C-12 


M 

Ti(3)  (k  ‘  1)2  Yi‘ 

k=0 


•M+k 


(C-13 


Also  from  the  definition 


( C  - 1 4 


We  can  now  use  Equation  C-7  to  calculate  the  non-recursive 
(truth)  sum  of  squared  residual  33R.  using  Equations  C-ll 
through  C-14. 


In  all  our  simulations  we  shall  calculate  the  following 
statistics  as  a  function  of  the  t-  (time  or  iteration  index 


(C-15) 


ABS  Error  1^1)  =  1^(1)  -  1^1)1 
ABS  Error  1.(2)  =  |1.(2)  -  7.(2) l  (C-16) 


ABS  Error  1i(3)  =  11,(3)  -  T.(3)J  (C-17) 


Percentage  error  SSR  = 


ABS  [SSRi  -  SSRi]  10Q 


(C-18) 


Figures  C-l  through  C-8  give  the  above  statistics  for  the  recur¬ 
sive  modification  using  single  precision  arithmetic  on  a  VAX  11/780 
computer.  The  amplitude  A  is  set  at  1.  The  period  T  is  set  at  100. 
other  parameters  are  varied  as  stated. 


The 


Note  for  all  simulations  and  for  all  parameter  settings,  the  modi 
fication  to  eliminate  the  accumulation  of  errors  in  the  recursive  LMSF 
procedure  described  in  Section  4.2  yields  identically  zero  errors  for 
all  statistics. 


Examining  Figures  C-l  through  C-8,  we  notice  a  moderate  growth  in 
absolute  error  of  the  linear  parameter  1(2)  and  a  substantially  larger 
growth  in  absolute  error  in  the  quadratic  parameter  1(3).  Overall, 
there  is  an  irratic  but  steady  growth  in  error  of  the  sum  of  squared 
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APPENDIX  D 


THEORETICAL  LOWER  RMS  ERROR  LIMIT 
FOR  THE  DYNAMIC  PROGRAMMING  ALGORITHM  DESCRETI2ED  POINTS 


From  the  above  diagram,  we  can  model  discretation  error  E  as  the 
following  linear  function 


E  *  t,  ■T<t<x  .  (0-1) 

This  model  is  repeated  many  times  over  many  discretized  points 
during  a  simulation. 


With  our  model  for  the  error  E  we  can  calculate  the  RMS  error  from 


the  definition  shown  below. 


The  above  evaluates  to  the  following. 

RMS  *  Vb2/12  .  (D- 

Now  for  our  simulations  B  =  0.25  (a  quarter  of  the  lag  spacing) 
Therefore,  the  limit  on  the  RMS  time  delay  difference  error  is  given 
below. 

RMS  =  0.0722  .  (D- 


-•  .*  / 
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